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Abstract 

This  report  summarizes  the  major  research  work  accomplished 
under  the  contract  including  the  topics  of  statistical  image  models, 
comparative  evaluation  of  image  processing  techniques,  image 
segmentation  algorithms,  two-dimensional  maximum  entropy  spectral 
analysis,  and  spatial  clustering  algorithms  with  applications  to 
artificial  and  remotely  sensed  images.  Detailed  list  of 
publications  available  in  open  literature  is  provided.  A  list 
of  software  package  generated  is  included  in  the  Appendices,  — 

I.  Introduction 

This  report  is  organized  according  to  the  topics  we  have  worked 
under  this  Contract.  A  brief  summary  of  each  is  presented.  Detailed 
list  of  publications  is  provided  in  References.  Over  40  technical 
reports  prepared  under  the  Contract  are  not  listed  however  as  most 
results  documented  in  reports  were  published  as  listed  in  References. 
Copies  of  two  papers  are  included  in  the  Appendices.  A  detailed 
list  of  image  processing  software  package  generated  is  also  included 
in  the  Appendices.  The  program  listings  in  magnetic  tapes  were 
delivered  to  Dr.  Doug  DeFriest  in  October  1981. 

II.  Major  Research  Results 

1.  Statistical  image  processing  techniques  for  additive  noise  case 
were  compared.  Median  filtering  followed  by  Kalman  adaptive 
filtering  is  most  effective.  For  Seasat  images  the  multicative 
noise  removal  is  considered  by  using  local  statistics.  (Refs.  16,20) 

2.  Statistical  image  segmentation  studies  include  the  extensive 
comparative  evaluation  of  supervised  and  unsupervised  segmentation 
techniques  for  texture  and  infrared  images.  The  segmentation  is 
performed  by  pixel  classification.  Both  Fisher's  linear  discriminant 
and  the  maximum  a  posteriori  estimation  procedures  are  found  to  be 
very  effective.  Statistical  techniques  however  are  limited  to  pixel 
based  segmentations.  (Refs.  4,9,10,12,13,14) 

3*  Statistical  image  modeling  study  is  concerned  with  the  auto¬ 
regressive  models  and  low  order  ARMA  models.  Such  modeling  leads  to 
image  enhancement,  segmentation  and  classification.  These  models 
provide  a  nice  way  to  take  into  account  the  contextual  dependence 
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among  the  nearest  neighbors.  The  question  remains  whether  the 
object  boundary  should  have  a  separate  model  from  the  remaining 
homogeneous  parts  of  the  image.  (Refs.  15 »18) 

4.  An  automatic  spatial  clustering  algorithm  has  been  developed 
for  image  segmentation  and  compression.  The  algorithm  can  determine 
the  minimum  number  of  clusters  and  can  also  work  with  a  specified 
number  of  clusters.  The  algorithm  has  been  successfully  tested  with 
various^  images  including  USC  image  data  base,  Seasat  images,  and 

17. S.  Army  topographic  images.  (Refs.  3*4) 

5.  A  two-dimensional  maximum  entropy  spectral  analysis  algorithm 
was  thoroughly  developed  and  tested  for  texture  image  analysis, 
classification,  segemntation  as  well  as  general  purpose  spectral 
computation  based  on  limited  number  of  data  points.  (Refs.  1,2,4,5» 

6,8) 

6.  An  initial  effort  of  tracking  image  sequence  was  made  by  using 
pixel  classification  for  object  extraction.  Further  study  to  model 
the  statistics  of  image  variation  is  much  needed.  (Ref.  4) 
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Abstract 

Effort  in  the  past  on  the  use  of  spectral 
features  for  discrimination  of  texture  images  has 
had  limited  success  and  the  feature  measures  com¬ 
puted  from  the  co-occurrence  matrix  are  preferred. 
In  this  paper  the  superiority  of  spectral  features . 
for  texture  classification  is  demonstrated.  A  new 
two-dimensional  maximum  entropy  spectral  analysis 
is  developed  which  provides  superior  resolution 
capability.  Thus  accurate  power  Bpectrum  can  he 
determined  from  which  various  ring  and  wedge 
spectral  features  are  computed  in  polar  coordinates. 
Extensive  conputer  results  reported  Indicate  that 
the  spectral  features  so  computed  provide  not  only 
a  good  measure  of  texture  coarseness  and  direction¬ 
ality,  but  also  conq>  arable  or  better  classification 
performance  than  that  reported  earlier.  A  typical 
performance  of  over  80  percent  correct  classifi¬ 
cation  is  available  from  the  extracted  spectral 
features  by  using  the  Fisher's  linear  discriminant 
for  classification.  A  set  of  nozmsLlized  features 
which  use  both  ring  and  wedge  features  is  particu¬ 
larly  recommended.  Computationally  the  method 
described  in  this  paper  is  far  better  than  the  use 
of  co-occurrence  matrix  as  the  iterative  algorithm 
used  for  spectrum  estimation  is  very  fast,  even 
with  the  use  of  a  minlconputer. 

I.  Introduction  * 


discrimination  with  the  energy  than  with  the  co¬ 
occurrence  measures . 

A  fundamental  problem  with  the  power  spectrum 
analysis  is  the  computational  accuracy  and  compu¬ 
tational  complexity.  For  texture  study,  accurate 
power  spectrum  must  be  computed  from  the  small 
image  segments.  In  this  case,  the  two-dimensional 
Fourier  analysis  cannot  provide  sufficient  accuracy 
as  the  Fourier  analysis  is  more  accurate  with  a 
large  number  of  pixels.  The  two-dimensional  maximum 
entropy  spectral  analysis,  however,  is  very  suitable 
for  a  small  number  of  pixels.  The  computational 
complexity  has  been  a  drawback  in  using  the  two- 
dimensional  maximum  entropy  spectral  estimation 
methods.  Recently,  Lim  and  Malik  [It]  have  proposed 
an  efficient  iterative  algorithm  for  the  two- 
dimensional  maximum  entropy  power  spectrum  esti¬ 
mation  which  can  obtain  good  resolution  and  suf¬ 
ficient  accuracy  for  the  finite  sample  two- 
dimensional  data.  A  study  of  the  Bpectral  esti¬ 
mation  of  texture  image  has  been  proved  to  be 
successful  [5]  by  using  a  minicomputer.  In  thiB 
paper,  this  method  is  used  for  the  calculation  of 
spectral  features  of  texture  image.  In  section  II, 
the  two-dimensional  maximum  entropy  power  spectrum 
estimation  is  briefly  discussed.  The  method  of 
selection  of  features  will  be  described  in  section 
III  while  section  IV  provides  some  experimental 
results  of  the  textural  classification. 


Although  it  is  generally  recognised  that 
texture  Images  contain  statistical,  spectral  and 
structural  domain  information,  the  use  of  spectral 
information  alone  can  be  quite  effective  in  the 
texture-image  analysis  studies  such  sb  texture  dis¬ 
crimination  and  segsietnation.  Bajcsy  and  Liberman 
[1]  expressed  the  power  spectrum  in  polar  coordi¬ 
nates,  then  integrate  over  r  and  4  to  obtain  the 
two-dimensional  functions.  The  location  of  peaks 
in  these  functions  indicates  prominent  texture 
coarseness  and  directionality.  Weszka  et.al.  [2] 
integrated  the  power  spectrum  within  16  spatial 
frequency  scjcs  which  were  combinations  of  four  1- 
octave  frequency  ranges  and  four  It 5°  orientation 
sectors.  They  also  computed  eight  "contrast" 
measures  based  on  the  cooccurrence  matrix,  and 
obtained  better  discrimination  than  with  the  power 
spectrum  measures .  Lows  [3]  computed  a  number  of 
energy  measures  by  filtering  the  texture  with  sets 
of  small  linear  operators,  then  squaring  and  summing 
the  output  of  each  filter.  He  reported  better 


Two-Dimensional  Maximum  Entror 
Spectrum  Estimation 


The  basic  concept  of  the  maximun  entropy  method 
(MEM)  of  spectral  estimation  is  to  extrapolate  the 
autocorrelation  function  of  a  random  process  by 
maximizing  the  entropy  H  of  the  corresponding  proba¬ 
bility  density  function 


H  ■  |  J  log  WgJdUjdiOg 


1  2”-w 


where  ^(u^,  u2)  is  the  power  spectrum  estimate  of 

the  random  process  x(n  ,  n^).  The  characteristics 
of  this  method  are  equivalent  to  the  autoregressive 
signal  modeling  and  the  power  spectrum  is  calculated 
by 
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P*(V  “2)  -  J  [X(ai,n2)e-J“lnl  e->2“2 
(nltn2)eA  (2) 

where  A(n^,  n^)  is  the  autocorrelation  vhose  power 
spectrum  is  l/Pfo^,  Ug)  and  A  is  a  set  of  points 
(n^,  n2)  where  the  autocorrelation  is  known. 

Since  the  filter  coefficients  cannot  be  ob¬ 
tained  directly  by  solving  the  normal  equation  as 
in  the  one- dimensional  case,  Lim  and  Malik  de¬ 
veloped  a  new  iterative  algorithm,  using  adaptive 
*  filtering  concepts.  The  basic  idea  of  this  algo¬ 

rithm  is  on  the  notion  that  the  given  correlation 
point  in  region  A  is  consistent  and  the  correspond- 
«  ing  coefficient  should  be  zero  outside  region  A  and 

proceed  this  iteration  repeatedly  until  an  optimal 
solution  is  obtained.  That  is,  for  a. given 
Rx(nl* ’n2>  tor  (ni*  determine  P  (m^  Wg) 

such  that  Px(u»- ,  «2)  satisfy  (2)  and 

Rx^nl*  n2*  *  F_1  [Px^“l*  “2^  for  ^ni*  n2^eA 

A  simple  flowchart  is  shown  in  Fig.  1.  We 
begin  with  some  Initial  estimate  of  X(n^,  n2), 

obtain  the  corresponding  correlation  function, 
correct  the  resulting  correlation  for  (n^,  UgjeA 

with  the  known  Rx(n^,  n2),  obtain  the  corresponding 
A(n^,  n2)  from  the  correct  corrleatlon  function, 
and  then  replace  the  restating  ng)  with  zero 

fOr  (n^,  n2)eA.  This  completes  one  iteration  and 
the  corrected  A(i;  j  ,  ,  ^)  is  a  new  estimate  of 
n2) 

In  Lim  and  Malik's  paper,  the  calculation  of 
the  autocorrelation  , ,  n2)  is  limited  to  the 

closed  analytic  form  especially  for  the  two- 
dimensional  sinusoids.  A  generalization  of  this 
method  and  the  application  to  a  two-dimensional 
real  data  have  been  discussed  by  Chen  and  Ioung[5). 
Even  for  a  small  number  of  missing  correlation 
points,  the  algorithm  can  still  provide  an  accurate 
spectrum..  Figure  2a  shows  the  spectrum  of  two 
sinusoids  (0.1234,  0.3 456),  (0.3125,  0.200)  in 
white  noise,  based  on  a  5x5  correlation  data  set, 
at  slgnl  to  noise  ratio  of  0.67.  With  the  corre¬ 
lation  points  (-1,  -1)  and  (l,  l)  missing.  Fig.  Zb 
shows  the  resulting  spectrum  which  is  nearly 
Identical  to  Fig.  2a. 

III.  Feature  Selection  and  Class fication  Method 

„  We  use  two  features  to  classify  the  texture 

images.  It  is  generally  recognized  that  a  coarse 
texture  will  have  a  high  value  of  power  spectrum 
near  the  origin  while  in  a  fine  texture,  the  value 
will  be  more  spread  out.  Thus,  if  one  wishes  to 
analyse  texture  coarseness,  a  set  of  features  that 
should  be  useful  are  the  averages  of  the  pover- 
speotnm  values  taken  over  a  ring-shaped  region 
centered  at  the  origin.  In  this  paper,  we  consider 
only  the  first  quadrant  of  the  power  spectrum,  then 


f*/2  . 
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for  various  values  of  r,  the  ring  radius. 

For  the  discrete  case,  this  can  be  vritten  as 
(for  rings  between  radius  r^  and  r^): 
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I  P»(x,y);  0  <.  x,y  <  n-1 

_2. _ 2 .  2,_2 


r  <x  +y  <r2 


(4) 


Similarly,  it  is  veil  known  also  that  the 
angular  distribution  of  pover  spectrum  is  sensitive 
to  the  directionality  of  the  texture  in  frequency 
m.  A  texture  with  many  edges  or  lines  in  a  given 
direction  8  will  have  high  values  of  power  spectrum 
around  the  perpendicular  of  8  +  */2;  while  in  a 
nondirectional  texture  the  spectrum  should  also  be 
nondlrectlonal.  Thus  a  good  set  of  features  for 
analyzing  the  texture  directionality  should  be  the 
averages  of  the  power-spectrum  values  taken  over 
wedge-shaped  regions  centered  at  the  origin,  i.e. 

♦e  “  f  P(r,  8)dr  (5) 

o 

for  the  various  values  of  8,  the  vedge  slope. 

For  the  discrete  case,  this  is  (the  vedge  be¬ 
tween  8^  and  8g)  given  by 

Vo"  Z.i  Vx’y)  (6) 

e^ftan  1  y/x<02 
0<x,y<n-l 


The  features  calculated  by  (4)  and  (6)  are  sensitive 
to  size  and  orientation  respectively,  but  not  to 
both.  In  order  to  obtain  the  comparable  feature 
sets,  we  obtain  a  set  of  equalized  features  by 
taking  the  average  over  the  intersection  area  of 
rings  and  wedges.  These  equalized  features  are 
also  studied  in  section  IV.  After  the  calculation 
of  features,  the  Fisher  discriminant  technique  is 
used  for  classification  [6] . 

IV.  Experimental  Results 

Because  of  the  computational  requirements  of 
the  method  and  the  limited  memory  capacity  of  the 
PDP  11/45,  all  test  samples  artL.  stored  in  our  DEC 
20  system  and  sent  through  a  communication  line  to 
the  PDP  11/45  for  the  spectrum  computation.  The 
test  samples  are  the  texture  images  taken  from  the 
USC  data  base.  To  verily  the  sensitivity  both  in 
coarseness  and  directionality,  we  select  some 
textures  that  contain  such  informations.  The  test 
samples  contain  six  classes  of  texture  (each  one 
has  four  samples)  and  are  shown  in  Fig.  3.  Each 
data  is  a  32x32  array  of  gray  level  0-255*  The 
pictures  of  class  1  reappear  buf  are  two  times 
larger  in  Fig.  4(a).  Figure  4(b)  is  the  correspond¬ 
ing  estimated  power-spectrum  display  of  the  upper 
left  data  in  each  class.  The  spectra  of  all  classes 
are  different  either  in  radial  or  angular 
distribution  [7]. 
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The  feature  sets  used  in  this  paper  are: 

*  r  for  (rl»  r2}  *  tl,3),  (3,6),  (6,12), 

1  2  (12,24) ,  (24^8) 

wedge:  for  (0^  0g)  -  (0,15),  (15,30), 

1  2  (30,45),  (45,60),  (60,75), 

(75,90) 

The  maximum  ring  radius  used  is  48  since  it  already 
covers  most  part  of  a  64x64  array  power  spectrum. 

A  combination  feature  of  ring  and  wedge  has 
been  tested  for  30  pairs  of  feature  values.  Table 
I  shows  part  of  features  which  did  higher  than  19 
out  of  24  correct,  i.e. ,  more  than  80*  correct. 
Table  II  shows  the  best  performing  pairs  using  the 
same  hind  of  features  (ring  and  ring,  wedge  and 
wedge):  there  are  6  out  of  25  pairs  which 
classified  correctly  higher  than  75*.  Other 
pairs'  results  are  concentrated  near  12-17,  i.e., 
more  than  50* -correct  recognition.  For  the  pairs 
that  contain  the  wedge  near  the  edges,  the  results 
are  very  good  Bince  the  test  samples  give  some 
directional  information.  Also  for  the  rings  a 
little  farther  from  the  origin,  the  results  are 
better  since  it  shows  a  large  difference  in  the 
spectrum  value  there. 

Equalized  features  are  also  tested:  we  used 
five  rings  intersected  with  three  wedges  (ring: 
(1,3),  (3,6),  (6,12),  (12,24),  (24,48)  and  wedge: 
(0,3fl),  (30, 60),  (60,90)).  105  pairs  of  features 

have  been  tested.  Table  III  shows  the  best  per¬ 
forming  pairs  of  which  the  best  score,  23  out  of 
24,  is  95*  correct.  From  the  results,  we  can  see 
that  the  ring  feature  (24,48)  gives  very  useful 
classification  information  indicating  that  there 
exists  a  large  textural  variation  in  that  region  aa 
the  texture  coarseness  plays  an  ioportant  role  in 
the  pair. 

V.  Discussion 

In  this  paper,  we  have  observed  that  equalised 
features  did  better  than  unequadized  ones  for  this 
test  samples.  It  is  shown  that  both  the  coarseness 
and  the  directionality  are  important  factors  in 
texture  discrimination.  For  the  consideration  of 
practical  use  in  automatic  classification,  various 
kinds  of  textures  must  be  tested  and  compared  with 
other  methods  using  the  non- spectral  features. 
Another  important  factor  which  may  influence  the 
results  is  that  if  we  increase  the  autocorrelation 
function  and  the  discrete  Fourier  transform  length 
while  estimating  the  power  spectrum,  the  accuracy 
and  the  resolution  will  be  better.  But  there  is  a 
tradeoff  between  the  accuracy  and  the  computational 
time.  In  this  paper,  these  parameters  (i.e.,  auto¬ 
correlation  function:  7x7,  discrete  Fourier  trans¬ 
form  length:  32)  are  chosen  for  the  real-time  pro¬ 
cessing  purpose.  Also  the  locations  of  the  main 
and  aecond  components  of  frequencies  can  serve  as 
another  important  features  because  they  vary  among 
different  textures. 
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Ring 

Wedge 

(21*  ,1*8) 

(0,15) 

22 

(1,3) 

(0,15) 

21 

(3,6) 

(0,15) 

20 

(1,3) 

(1*5,60) 

19 

(12,21*) 

(0,15) 

19 

(3,6) 

(30,1*5) 

19 

(3,6) 

(75,90) 

■>o 

Table  I:  Best  perfoming  pairs  using  the 

combination  feature  of  ring  and  wedge 
for  those  with  more  than  80%  correct 
classification. 


Features 


(6,12) 

(6,12) 


(30,45) 

(15,30) 

(30,1*5) 

(1*5,60) 


Number  correctly 
classified 


Ring 

(21*.  1*8) 

20 

(12,21*) 

19 

Wedge 

(75,90) 

20 

(60,75) 

18 

(60,75) 

18 

(60,75) 

18 

Table  II:  Beat  performing  pairs  using  same  kind 
of  features,  for  those  with  more  than 
75 t  correct  classification. 


Features  Number  correct 

classification 


Bing  n 

Wedge 

Ring  n 

Wedge 

(3,6) 

(0,30) 

(2l*,l*8) 

(0,30) 

23 

(12, 21*) 

(60,90) 

(2U.U8) 

(0,30) 

23 

(2U.U8) 

(0,30) 

(2l*,l*8) 

(30,60) 

22 

(1,3) 

(0,30) 

(2l*,l*8) 

(0,30) 

22 

(1,3) 

(30,60) 

(2l*,U8) 

(0,30) 

22 

(1,3) 

(60,90) 

(24, 1.8) 

(0,30) 

22 

(3,6) 

(30,60) 

(24,48) 

(0,30) 

21 

(3,6) 

(60,90) 

(24,48) 

(0,30) 

21 

(6,12) 

(0,30) 

(24,48) 

(0,30) 

21 

(12,21*) 

(30,60) 

(24,48)  - 

“(0,30) 

21 

(12,21*) 

(0,30) 

(12,24) 

(60,90) 

20 

(1,3) 

(30,60) 

(6,12) 

(0,30) 

20 

(1,3) 

(60,90) 

(6,12) 

(0,30) 

20 

(12,21*) 

(0,30) 

(24,48) 

(0,30) 

20 

(1,3) 

(0,30) 

(6,12) 

(0,30) 

19 

(3,6) 

(0,30) 

(6,12) 

(0,30) 

19 

(3,6) 

(30, 60) 

(6,12) 

(0,30) 

19 

(6,12) 

(0,30) 

(12,24) 

(60,90) 

19 

Table  III:  Best  perfoming  pairs  using  equalized 
features  for  those  with  more  than  80S 
correct  classification. 
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Appendix  C 


A  computationally  efficient  spatial  clustering  algorithm  is  presented  for  image 
segmentation  and  compression.  The  algorithm  can  automatically  determine  the  minimum 
number  of  clusters  and  can  also  work  on  a  specified  number  of  clusters.  Examples  are 
given  on  the  processing  of  Seasat  images  using  the  algorithm. 


Introduction 


Besides  the  use  of  acoustic  sensors,  remotely  sensed  images  can  provide  essential 
information  on  the  object  extraction  and  tracking  in  the  ocean  environment.  Seasat 
SAR  images  are  a  good  example.  Many  automatic  image  analysis  algorithms  have  been 
developed.  Such  algorithms  are  generally  application  dependent.  For  remote  sensing 
images,  cluster  analysis  is  important  as  it  reveals  the  structure  of  the  data  from 
which  useful  information  can  be  derived.  Conventional  clustering  methods  do  not 
preserve  the  spatial  relations  in  a  image.  Spatial  clustering  for  image  analysis  has 
been  considered  [l][2j.  However  feature  extraction  was  not  taken  into  account. 
Furthermore  the  computation  involved  is  quite  extensive.  A  more  efficient  spatial 
clustering  algorithm  is  developed  for  minicomputer  processing,  that  employs  properly 
selected  features.  The  clustered  image  shows  various  regions  (segments)  from  which 
desired  objects  may  be  extracted.  Furthermore  considerable  image  data  compression 
is  accomplished  with  essentially  no  loss  of  information.  Examples  are  based  exclu¬ 
sively  on  Seasat  images  dealing  with  the  ocean  environment. 


Algorithm  for  Spatial  Clustering 

The  algorithm  proceeds  as  follows: 

(1)  Form  a  feature  set,  for  each  pixel,  consisting  of  local  mean  and  gradient. 

Other  features  may  also  be  used. 

(2)  For  each  2x2  subarea,  measure  the  mean  vector  and  dispersion. 

(3)  Determine  the  critical  dispersion,  and  calculate  the  merging  distance  d. 

(4)  Merge  adjacent  subareas  with  distance  less  than  d  to  form  subregions.  Calculate 
the  mean  vectors  of  subregions. 

(5)  Group  these  mean  vectors  into  clusters  using  K-mean  algorithm  which  converges 
to  several  cluster  centers  representing  the  mean  vectors  of  regions. 

For  a  given  inter-region  threshold  distance,  the  algorithm  can  automatically 
adjust  to  an  appropriate  number  of  clusters.  If  the  number  of  clusters  is  specified 
as  in  conventional  cluster  analysis,  the  algorithm  will  provide  the  specified  number 
of  clusters. 


The  image  speckles  in  synthetic  aperture  are  multiplicative  in  nature.  Without 
removing  such  noise,  the  cluster  results  may  still  be  very  noisy.  A  simple  pre¬ 
processing  method  is  to  use  the  Sigma  filter  suggested  by  J.S.  Lee  [3],  For  each 
5x5  (or  3x3)  subarea  the  average  gray  level  value  is  compared  with  the  three-standard 
deviation  (3o)  of  the  normalized  image  histogram  (first  order  probability  density) . 

If  the  value  is  within  3o  tLen  the  center  pixel  is  replaced  by  the  average  value.  If 
the  value  exceeds  3a,  then  it  is  an  indication  of  edges  or  object  boundary  and  the 
original  gray  level  is  retained.  The  procedure  thus  provides  a  compromise  between 
noise  filtering  and  edge  preserving  and  adds  only  slight  amount  of  computation  to 
the  clustering  algorithm. 

Computer  Results 

The  original  Seasat  images  are  all  of  256  gray  levels.  For  convenience  with 
minicomputer  processing  the  digitized  pictures  are  all  reduced  to  256x256  pixels 
even  though  the  original  images  are  much  larger  in  size.  The  cluster  algorithm  is 
set  such  that  a  maximum  of  7  clusters  is  selected.  Figure  1  is  for  the  scene  of  a 
ship  off  Chesapeake  Bay.  Figure  la  is  the  original.  Figure  lb  is  the  spatially 
clustered  image.  Figure  lc  has  the  histograms  of  original  (left)  and  clustered 
(right)  images.  Figure  Id  uses  the  Sigma  filter  preprocessing  with  3x3  subarea 
while  Fig.  le  is  the  corresponding  result  with  5x5  subarea.  Preprocessing  with  5x5 
subarea  followed  by  spatial  clustering  appears  to  be  the  best.  Figure2  is  the 
scene  of  Anchorage,  Alaska  with  Fig.  2a  for  original  and  Fig.  2b  for  5x5  preprocess¬ 
ing  followed  by  spatial  clustering.  Figure  3  is  the  scene  of  Nantucket  Shoals  with 
Fig.  3a  for  original,  Fig.  3b  for  5x5  median  filtering  and  Fig. 3c  for  5x5  pre¬ 
processing  followed  by  spatial  clustering.  Computer  results  all  show  that  the 
"natural"  clusters  of  the  original  images  are  very  much  preserved  while  noises  are 
considerably  reduced,  and  at  the  same  time  the  contrast  is  enhanced,  with  the  use  of 
the  spatial  clustering  algorithm. 
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Appendix  D 


APPENDIX  PRESPA. FOR 


c  ****************************************** 

C  PRESPA. POR  <  HAS  PREA4.F0R  ) 

C  30-MAY -82  (  300  SECTION  :  READ  DATA  FILE  FROM  TAPE  ) 

C  27-APR-82 

C  PROGRAM  TO  READ  A  DISK  FILE  256X256  PIXELS 

C  OUTPUT  ITS  FEATURE  VECTORS,  MEAN  ,  GRADIENT 

C  OR  ORIGINAL  INTENSITY  ,  GRADIENT 

C  CHOICE  3:  READ  DATA  FROM  TAPE 

C  TO  FORM  NOF1  ,  BOP2  OUTPUT  FILES 

C  FEATURE  1  AND  FEATURE  2  RESPECTIVELY 

C  NOFl:  OUTPUT  DATA  ,  MEAN  OR  ORIGINAL  GRAY  LEVEL,  COMPONENT 
1 

C  NOF2:  OUTPUT  DATA  ,  GRADIENT,  COMPONENT  2 

C  IF  CHOICE  1  OR  2  S  READ  DISK  FILE  MOF4 

C  NOF4:  INPUT  DATA 

C  ******************************************* 

INTEGER  F (1024), F2 (102 4), CHOICE, IMEAN(256),IGRAD(256) 

REAL  DMEAN(256),DGRAD(256 ), MSI (256),HS2(256) 

DATA  NOFl,NOF2, NOP 4/1, 2,4/ 

1001  FORMAT! *  PROGRAM  PRESPA.FOR*/ 

2  *  PREPROCESSING  IMAGE  DATA  TO  FORM  FEATURE 
VECTORS' / 

3  *  FOR  AUTO  SPATIAL  CLUSTERING'/ 

4  *  INPUT:  FTN4.DAT  OR  TAPE  -  ORIGINAL  GRAY  LEVEL*/ 

5  '  OUTPUT:  FTN1.0AT  FEATURE  COMPONENT  1*/ 

6  '  FTN2.DAT,  FEATURE  COMPONENT  2*/ 

7  '  choices:*/'  i:  local  mean  ,  local  gradient'/ 

8  *  2:  ORIGINAL  INTENSITY  ,  LOCAL  GRADIENT*/ 

9  *  3:  COMPRESS  A  2 BY  2  SUB  IMAGE  INTO  1  PIXEL*/ 

1  *  4:  READ  TAPE  DATA  (NORGL*NORCP)  TO 

DISK(NOL*NOP )*/) 

1002  FORMAT! 15) 

1003  FORMAT!*  ENTER  INPUT  AND  OUTPUT  FILES  SIZE*/ 

1  *  NOLIN,NOPIN,NOL,NOP:  FORMAT! 415) *) 

1004  FORMAT  (415) 

1  CONTINUE 

HRITE!7,1001) 

R E AD(5, 1002 )CHO ICE 

IF  (CHOICE. LE.O. OR. CHOICE.GT.4)  GOTO  1 
GOTO  (50,50,300,400), CHOICE 
50  CONTINUE 

HRITE(7,1003) 

READ(5,1004)NOL IN, NOPIN, NOL, NOP 
C  ORIGINAL  INTEGER  DATA  FILE  OF  AN  IMAGE 

DEFINE  FILE  NOF4(NOLIN,NOPIN,U,INDX4) 

NOLl=NOL-l 

NOPUNOP-1 

C  OUTPUT  FILES,  NOFl  ,  NOF2,  INTEGER  NUMBERS 

DEFINE  FILE  NOFl (NOL, NOP, U, INDX1) 

DEFINE  FILE  NOF2(NOL,NOP,U, INDX2) 

GOTO  (100, 200), CHOICE 
100  CONTINUE 
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70 


90 


200 

C 

C 


270 


290 


300 

C 

c 

c 

c 

1031 

TAPE* 


1032 

1033 


1034 


DO  90  Isl,NOLl 
INDX4-I 

R£AD(M0F4*INDX4)(F(K),K=1,N0PIN) 

RE AD(N0F4*IMDX4)(F2(K)/K-lsN0PlN) 

DO  70  J=1,N0P1 

IMEAN(J)=(F(J)*F(J+l)+F2(J)+F2(J+l))/4 

IGRAD(J)=(IABS(F(J)-F(J+l))+IABS(F2(J)-F(J)))/2 

CONTINUE 

IMEAN(NOP)=IMEAN(NOPl) 

IGRAD( HOP )~IGRAD(N0P1 ) 

INDXl-I 

INDX2-I 

HR1TE( NOFl'INDXDOMEAN (K)*K=1,NQP) 
NRITE(NOF2*INDX2)(IGRAD(K),K=l,NOP) 

CONTINUE 

NRITE(NOFl'INDXl)(IMEAN(K),Ksl,NOP) 

WRITE(NOF2*INDX2)(IGRAD(K)sK=l,NOP) 

GOTO  900 
CONTINUE 

CHOICE  2S  SECTION 

FEATURE:  ORIGINAL  INTENSITY  ,  LOCAL  GRADIENT 

DO  290  1=1, NOLI 

INDX4-I 

READ(N0F4*INDX4)(F(K),K=1#N0P IN) 

RE AD(N0F4*IMDX4)(F2(K), K=l, NOPIN) 

DO  270  J=l,NOPl 

IGRAD(J)=(IABS(F(J)-F(J+l))+IABS(F2(J)-F(J)))/2 

CONTINUE 

IGRADC  NOP )=IGRAD(N0P1 ) 

INDX1=I 

INDX2-I 

WRITE! N0F1*INDX1)(F(K),K=1, NOP) 
WRITE(NOF2*INDX2)(IGRAD(K),K=l,NOP) 

CONTINUE 

NRITE(N0F1*INDX1)(F2(K),K=1,N0P) 

WRITE! NOF2*INDX2)( IGRAD(K),K=1, NOP) 

GOTO  900 
CONTINUE 

READ  TAPE  DATA  FILE  SECTION 

NOLI#NOLF#NOPI,NOPF S  COVER  THE  AREA  INTERESTED 
NOL,  NOP:  SIZE  OF  THE  OUTPUT  DATA  IN  FEATURE  SPACE 

EACH  RECORD  REPRESENTS  2  BY  2  ORIGINAL  PIXELS 
WRITE(7,1031) 

FORMAT! *  ENTER  ORIGINAL  TAPE  FILE  SIZE  AND  WHICH  FILE  IN 

/ 

1  *  NORGL,NORGP,NTH  (  FORMAT(3I6)  ):*> 

READ(5,1032)NORGL,NORGP,NTH 

FORMAT (31 6) 

NR ITEC 7,1033) 

FORMAT ( *  WHICH  PART  OF  IMAGE  TO  BE  PROCESSED?'/ 

1  *  NOLI,NOLF,NOPI,HOPF:  (  F0RMAT<4I6>  )*) 

READ(5#1034)NOL I,NQLF,NOPI,NOPF 
FORMAT (416) 

NOLs(NOLF-NOLI*l)/2 

N0Ps(N0PF-N0Pl4l)/2 
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UR1TE(7,1035)NOLI,NOLF,NOPI,NOPF,NOL,NOP 
1035  FORMAT!*  CHECK:  NOLI,NOLF,NOPI,NOPF, MOL, NOP'/ 

1  616 

2  /*  NOL  X  NOP  HILL  BE  THE  OUTPUT  SIZE*) 

DEFINE  FILE  NOFl(NOL,NOP,U, INDX1 ) 

DEFINE  FILE  NOF2(NOL,MOP,U,IMDX2) 

REHIND  NOF1 
REHIND  NOF2 

C  SKIP  THE  PART  NOT  TO  BE  PROCESSED 
NSKIP=NOLI-l 

C  NOLI  IS  THE  FIRST  LINE  TO  BE  PROCESSED 
IF  (NSKIP.LT. 1)  60T0  312 
DO  310  I=1,NSKIP 
310  CALL  READUN(F,1,NTH) 

312  CONTINUE 

DO  350  1=1, NOL 
DO  315  J=l,NOP 
IMEAN(J)=0 
IGRAD(J)=0 
315  CONTINUE 

CALL  READUN(F,1,NTH) 

CALL  READUN!F2,1,NTH) 

DO  320  J=l,NOP 

IMEAN(J)=(F(N0PI*J*J-2>4F(N0PI>Kl-frJ-l)*F2(N0PI-»J^J-2)* 

1  F2(NOPI*J*J-l))/4 

I GRAD( J )=(IABS(F( NOPI+J+J-2 )-F (NOPI* J* J-l ) )♦ 

1  IABS(F(N0PI+j4j-2)-F2!N0PI+J+J-2)))/2 

320  CONTINUE 
INOX1-I 
INDX2=I 

WRITE! M0F1*INDX1)( 1MEAN !K),K=l,NOP) 

WRITE! NOF2'INDX2)(IGRAD(K),K=l, NOP) 

350  CONTINUE 
GOTO  900 
400  CONTINUE 

C  READ  TAPE  DATA  FILE  TO  FORM  A  COMPRESSED  DISK  DATA  FILE 
C  INPUT  SIZE  :  NORGL  BY  NORGP  PIXELS 
C  OUTPUT  SIZE:  NOL  BY  NOP 
URITE(7,1041) 

1041  FORMAT (  *  ENTER  NORGL,NORGP,NOL,NOP,NOF,NTH*/ 

1  *  E.G.,  256,256/64/ 64,?/?*) 

RE  ADC 5 , 10  42) NORGL, NORGP, NOL, MOP, NOP, NTH 

1042  FORMAT (616) 

DEFINE  FILE  NOF(NOL,NOP,U,INDEX) 

IT ERL- NORGL/ NOL 
IT ERP- NORGP/ NOP 
DO  450  1=1, NOL 
DO  420  J=l, ITERL 
420  CALL  READUN(F,1,NTH) 

WRITE! NOF'INDEX )(F( ITER P*K),K=1, NOP) 

450  CONTINUE 
GOTO  900 
900  CONTINUE 
CALL  EXIT 
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C  9A123456789B123456789C123456789D1234567 89E123456789F12345678 

9G12 

C  FILE  NAIffi:  SPABH.FOR 

C  AUTO  SPATIAL  CLUSTERING  FOR  BLACK/WHITE  IMAGE  DATA 
C 

C  REFERENCE  FUKADA,"  SPATIAL  CLUSTERING  PROCEDURESFOR  REGION 

C  ANALYSIS1*, PATTERN  RECOGMITION,12,395-403,(1980). 

C 

C 

C  INPUTS:  NOFl,  FTN1.DAT,  FEATURE  1 

C  MOF2,  FTN2.DAT,  FEATURE  2 

C  OUTPUTS :N0FC,  FTN12.DAT,  CLUSTERING  RESULT  FOR  COLOR 

DISPLAY 

C  NOFB,  FTN11.DAT,  CLUSTERING  RESULT  FOR  BLACK/WHITE 

C  DISPLAY 

C  NOFF,  FTN15.DAT,  SOME  PARAMETERS  DURONG  PROCESSING 

C  SIZE  OF  IMAGE  IS  RESTICTED  TO  256  BY  256  IN  FEATURE  SPACE 
C  - 

LOGICAL*!  DUMMY (9) 

REAL  PREV(60,2),PRENO(60) 

COMMON  /BLQCKO/IOP (10), NOFl, NOF2,NOFB,MOFC, NOFF 
COMMON  /BLOCK 1/NOL, NOP, NO L2,NOP2,NOLH, NOP H 
COMMON  /BLOCK2/CKV(60,2),CNO(60),FKV (30,2),FNO(30) 

COMMON  /BLOCK3/IA(256),IB(256),A(256,2),B( 256,2) 

COMMON  /BL0CK4/ AB (256,2 ),ABM( 2, 128,2 ), TRACE (128) 

COMMON  /BLQCK5/LABEL(2,l28),ICQLOR(256),IBU(256) 

OATA  NOFl, NOF2, NOFB, NOFC,NOFF/l,2,ll,12,15/ 

1701  FORMAT ( *  RUNNING  PROGRAM  SPABH.FOR  AUTO  SPATIAL 
CLUSTERING*/ 

1  *  INPUTS  ,  OUTPUTS:  NOFl,NOF2,NOFB,NOFC,NOFF'/ 

6  *  FOR  CHECK: *,515/ 

2  *  NOFC  ,  NOFB  STORE  TEMPORARY  DATA  DURING 
PROCESSING*/ 

3  *  OUTPUTS:  NOFC (REWRITTEN)  7  OR  LESS  COLORS*/ 

4  '  NOFB (REWRITTEN)  BLACK  AND  WHITE 

DISPLAY*/ 

5  *  NOFF,  INFORMATIONS  DURING  PROCESSING*) 

1515  FORMAT ( *  MERGE  ITERATION*, 15 ) 

1520  FORMAT!*  THE  *,I5,*-TH  ITERATION  REACHES  MAXIMUM  NO. 
CLUSTERS*/ 

1  *  NO.  OF  CLUSTERS:*, 15) 

1511  FORMAT! *  KERNEL  CANDIDATE  VECTORS  FOR*, 14,*  STARTING  *, 

1  'CLUSTER  CENTERS*/) 

1512  FORMAT!*  ENTER  IMAGE  DATA  FILE  SIZE  (  FEATURE  SPACE  )*/ 

1  *  NOL,  NOP:  F0RMAT(2I 5) *) 

1513  FORMAT(2I5) 

1514  FORMAT!*  CHECK: NOL, NOP, NOL2, NOP2,NOLH,NOPH*/6I5) 

1516  FORMAT ( *  ENTER  OPTIONS:  (IQP(K),K=1,10)*/ 

1  *  IOP(l):  CONTROLS  PRINTER,  1:  MEANS  ,  TRACE 
MATRIX'/ 

2  '  IOP(2):  K-MEANS  ALGORITHM  DETAILS  ON  SCREEN*/ 

3  *  IOP(3):  MERGE  DETAILS,  LABEL(2,K)  ARRAY*/ 


o  o  o  o  n 
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4  *  IOP<4):  NSTEP,  MO.  OP  STEPS*/ 

5  *  IOP<5):  1,  SKIP  K-MEAN  ITERATION'/) 

1517  FORMAT (10 16) 

1518  FORMAT! *  TODAY  IS  *,9A1) 

REWIND  HOFF 

WR  IT  E(  7,1501) 

WRITE! HOFF# 1501) 

CALL  DATE(DDMMY) 

WRITE! NQFF,1518)(DDMMY(K),K=1,9) 

WR IT E( 7, 151 8)  (DDMMY(K),K=1,9) 

WRITE! NOFF, 1701 )N0F1,N0F2, NOFB, NOFC, NOFF 
WRITE! 7,1701 )NOPl,NOF2, NOFB,NOFC,NOFF 
WRITE!7,1512) 

WRITE! NOFF, 1512) 

READ!5,1513)NOL,NOP 
WRITE!  7,1516) 

WRITE! NOFF, 1516) 

READ!5,1517)!IOP!K),K-1,10) 

WRITE! NOFF, 1517 )(IOP(K),K=l, 10) 

WR ITE( 7,1517) ( IOP(K),K=l, 10 ) 

NOL2=MOL+WOL 

NOP2=NOP+NOP 

NOLH-NOL/2 

NOPH-HOP/2 

WRITE! NOFF, 1514 )NOL, NOP, NOL2,NOP2,NOLH, NOPH 
WRITE! 7, 1514) NOL, NOP, NOL2,NOP2,NOLH, NOPH 
DEFINE  FILE  NOFl!NOL,NOP,U,INDXl) 

DEFINE  FILE  NOF2!NOL,NOP,U,INDX2) 

DEFINE  FILE  NOFB!NOL,NOP,U,INDXB) 

DEFINE  FILE  NOFC!NOL,NOP,D,INDXC) 

C  NOFF-FTN15.DAT  UNFORMATTED 
REWIND  NOF1 
REWIND  NOF2 
REWIND  NOFB 
REWIND  NOFC 

1501  FORMAT!*  THIS  IS  THE  LOG  FILE  OF  EXECUTING  SPABW.FOR*) 

C  - 

CALCULATE  MEANS  OF  FEATURE  VECTORS  OF  2  BY  2  SUBIMAGE 
STORE  IN  NOFB:  FIRST  HALF  —  FEATURE  1  MEANS  OF  128X126 

SECOND  HALF  —  FEATURE  2  MEANS 
128  X  128  SUBIMAGES  EACH 

IN  NOFC:  FIRST  AND  SECOND  HALF  ARE  THE  SANE,  TRACES 
CALL  DISPER 

C  FIND  MAX  ,  MIN  OF  TRACE  MATRIX 
CALL  MAXMIN(DMAX,DMIN) 

C  MERGING  SECTION 

NSTEPsIOP(4) 

STEP=(DMAX-DM  ID/FLOAT  (NSTEP) 

WRITE! 7,1522 )DMAX,DM1N, STEP 
WR ITE! NOFF, 1522 )DMAX,DMIN,STEP 
1522  FORMAT!*  DMAX**,E20.8, *  DMIN=*,E20.8,*  STEP**,E20. 8) 
IPREV-0 

C  ITERATIONS  TO  FIND  MAXIMUM  NO.  OF  CLUSTERS 
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DO  300  Isl,NSTEP 
IM=I 

URITE!NOFF,1515)IM 

CALL  MERGE! IM,STEP,DMIN,NCLSR) 

C  ACCEPTED  10.  OF  CLUSTERS:  7  OR  LESS 

IF  !IPREV.LE.7. AND.IPREV.GT. 1)  GOTO  333 
IPREV-NCLSR 

C  SAVE  CURRENT  NUMBER  OF  CLUSTERS  AND  KERNEL  VECTORS 
DO  200  Jsl,IPREV 
PRENO(J)=CNO(J) 

PREV!J,1)=CKV!J,1) 

PREV!J,2)sCKV!J,2) 

200  CONTINUE 
300  CONTINUE 
333  NI=IM-1 

WRITE! HOFF, 1520 )NI,IPREY 
DO  350  Jsl,IPREV 
CNO(J)=PRENO(J) 

CKV(J,1)=PR£V(J,1) 

CKV(J,2)=PREV(J,2) 

350  CONTINUE 

WRITE! N OFF, 1561 ) 

1561  FORMAT! *  MERGE  ENDED  WITH  MAXIMUM  NO.  CLUSTERS') 
WRITE!NOFF,1562)!!CKV!N/L),L=l,2),N=l,IPREV) 

1562  FORMAT!'  BEFORE  SORTING'/ 

1  '  KERNEL  CANDIDATE  VECTORS'/ 

2  30((5X,2E20.8)/)) 

C  SORT  THE  CANDIDATE  VECTORS 

NC=IPREV 

C  SORT  THE  CANDIDATE  KERNEL  VECTORS 

CALL  SORT!NC) 

NRITE!NOFF,1563)!!CKV!N,L),L=l,2),N=l,NC) 

1563  FORMAT!*  SORTED  KERNEL  CANDIDATE  VECTORS*/ 

1  30((5X,2E20.8)/)) 

IP  !NC.GT.7)  NC=7 

C  FOR  THE  PURPOSE  OF  AED-512  PSEUDO  COLOR  DISPLAY 
WRITE!NOFF,1511)NC 

C  IF  !I0P!5).EQ.l)  SKIP  THE  K-MEAN  ITERATIONS 
C  DIRECTLY  USE  MERGING  RESULT  CAD  ID ATE  KERNEL  VECTORS 
C  TO  CLASSIFY  THE  IMAGE 

WRITE!7,1568)IOP!5) 

1568  FORMAT!*  IOP!5)=',I5) 

IF  !I0P!5).NE.l)  GOTO  700 
WR1TE!7,1570) 

1570  FORMAT!*  SKIP  K-MEAN  ITERATION') 

DO  650  NS1,NC 
DO  640  L-1,2 
FK V ! N, L )*CKV !  N,L) 

640  CONTINUE 

650  CONTINUE 

GOTO  800  . 

700  CONTINUE 

WRITE!  7,1 580) 

1580  FORMAT!*  CALLING  KERVEC:  K-MEAN  ITERATION') 
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C  ITERATIONS  TO  FIND  MORE  ACCURATE  KERNEL  VECTORS 

C  CALLED  FINAL  KERNEL  VECTORS 

CALL  KERVEC(NC,KK,DD) 

HR1TEC  HOFF, 1500  >KK, DD 

1500  FORMAT! IX, 'CLUSTERING  REPEATS', 1X,I3,1X,'TIHES '/IX, 
l'TBE  FINAL  VITHIR-CLASS  DISTANCE  IS',lX,E20.8/> 

800  CONTINUE 

C  CLASSIFICATION  SECTION 

C  OUTPUTS:  NOFC,  COLOR  DISPLAY  RESULT 

C  NOFD,  BLACK/HHITE  DISPLAY  RESULT 

CALL  CL  ASS (NC) 

MRITE(M0FF,1523) 

1523  FORMAT (10 X, ' 111  COMPLETE  EXECUTION  OF  PROGRAM  SPABH  III') 
CALL  EXIT 
END 
C 
C 

C  -  SUBPROGRAMS  - 

C 

C  - 

C  SUBROUTINE  TO  CALCULATE  TRACE  MATRICES  OF  FEATURE  MATRICES 
C  STORED  IN  NOFD 

SUBROUTINE  DISPER 

COMMON  /BLOCKO/ IOP ( 1 0 ) , NOF1 , N OF2 , NOF B, NOFC, MOFF 
COMMON  /BLOCKl /MOL , NOP , NOL2 , NOP 2, NOLH, NOPH 
COMMON  /BLOCKS/ IA(256),XB(256),A( 256,2 ),B(256,2) 

COMMON  /BLOCK4/AB( 256,2 ),ABH(2, 128,2 ),TRACE(128) 

REWIND  NOF1 
REMIND  NOF 2 
REMIND  NOFB 
REMIND  NOFC 

C  PROCESS  THROUGH  ROMS  OF  DATA  MATRIX 

DO  100  Isl,NOLH 
12=1+1-1 
III  0X1=  1 2 
INDX2=I2 
DO  40  JJ=1,2 

C  READ  2  LINES  OF  EACH  FILE 

R£AD(N0F1*INDX1 )(IA(K),K=l,NOP) 

DO  10  Jsl,NOP 

10  A(J,4J)zFL0AT(IA(J)) 

READ(N0F2*INDX2)(IB(K),K=1,N0P) 

DO  20  «J si, MOP 

20  B(J,JJ)sFLOAT(IB(J)) 

40  CONTINUE 

C  CALCULATION  THROUGH  EACH  SUBIMAGE 

DO  80  K*l,NOPH 
K1=K+K-1 
K2*K1+1 
S1*0. 

32*0  m 

DO  62  M»K1,K2 
DO  60  L*l,2 
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S1=S1*A<M,L) 

S2-S2*B(M,L) 

60  CONTINUE 

62  CONTINUE 

ABN(2,K,1)=S1*0.25 
ABM(2/K/2 )=S2*0 .25 
S1=0. 

S2=0. 

DO  72  M-K1,K2 
DO  70  L=l,2 

S1=S1+!A!M,L)-ABM!2,K,  1))**2 
S2=S2*!B!M,L)-ABM!2,K,2))**2 
70  CONTINUE 

72  CONTINUE 

TRAC£!K )-!Sl+S2  )*0. 25 
80  CONTINUE 

INDXB>I 

IF  (IOP(l).EU.l)  WRITE(6,1001)(ABM(2,K,1),K=1,32) 
WRITE! NQFB'IKDXB) !ABM!2#Ksl )*K=ls NOPH) 
IMDXB=I+NOLH 

IF  (lOP(l).EQ.l)  WRITE!6,1002)!ABM!2,K,2),K=1,32) 
BRITE!NOFB*INDXB)!ABM!2,K*2),K=l,NOPH) 

INDXC=I 

WRITE! NOFC*INDXC)!TRACE!K),K=l,NOPH) 

INDXC= I+NOLH 

IF  (IOP(l).EQ.l)  WRITE!6,1004)!TRACE!K),K=1,32) 
WRITE! NQFC'INDXC) (TRACE (K)/K=l^ NOPH) 

100  CONTINUE 

1001  FORMAT! *  AM**32F4.0) 

1002  FORMAT!*  BM*,32F4.0) 

1004  FORMAT!*  TR*,32F4.0) 

RETURN 

END 

C  READ  TRACE  MATRIX  TO  FIND  MAX  ,  MIN 
C 

SUBROUTINE  MAXMIN!DMAX# DMIN) 

COMMON  /BLOCKO/IOPUO),NOF1,NOF2,NOFB,NOFC,NOFF 
COMMON  /BLOCKl/NOL,NOP,NOL2,NOP2,NOLH,NOPH 
COMMON  /BLOCK3/IA!256),IB!256),A!256,2),B!256,2) 
COMMON  /BLOCK4/AB!256,2),ABM!2/128,2),TRACE!128) 
REMIND  NOFC 
INDXC-1 

READ!N0FC*1NDXC)! TRACE (K)/K=l/ NOPH) 

DMAX=TRACE(1 ) 

DMIN=TRACE(1) 

DO  10  Js2,NOPH 

IF  !TRACE\J).LT.DMIN)  DMIN=TRACE! J) 

IF  !TRACE!J).GT.DMAX)  DMAX=TR ACE! J) 

10  CONTINUE 

DO  100  I=2/NOLH 
INDXC-I 

READ!NOFC*INDXC )! TRACE !K)/K=1/ NOPH) 

DO  30  J=l,NOPH 

IF  !TRACE!J).LT.DMIN)  DMIN=TRACE! J) 
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IF  (TRACE(J).GT.DMAX)  DMAX=TR ACE( J) 

30  CONTINUE 
100  CONTINUE 

WR IT E! 7,1 00 1 ) DM AX, DM IN 

1001  FORMAT!/'  DMAX=',F12.4,'  DMIN= *,F12.4) 

RETURN 

END 

C 

C  - 

C  MERGE  AND  DECIDE  KERNEL  CANDIDATE  VECTORS 

C 

SUBROUTINE  MERGE! IMRGE,DSTEP,DMIN,LLBS) 

REAL  F 1 RST! 60) 

COMMON  /BLOCKO/IOP!10),NOF1,NOF2,NOFB,NOPC,NOFF 
COMMON  /BLOCKl/NOL,NOP, NOL2,NOP2,NQLH,NOPH 
COMMON  /BLQCK2/CKV(60,2),CNO!60),FKV!30,2),FNO!30) 
COMMON  /BLOCK3/IA!256),IB!256),A!256,2),B!256,2) 
COMMON  /BLOCK4/AB!256,2),ABM(2,128,2),TRACE(128) 
COMMON  /BLOCK5/LABEL!2,128),ICOLOR!256),IBM!256) 

1502  FORMAT! 30!/*  CNO! N)  : ',F12.1  )  ) 

1503  FORMAT!'  *FMIN,  DMRGE: ',2E20 . 8 ) 

1504  FORMAT!/'  IMRGE:',I5) 

1505  FORMAT!'  FIRST  SUBIMAGE:') 

1506  FORMAT! '  JOIN  CLUSTER  NL:') 

1507  FORMAT!'  NEN  CLUSTER:') 

1508  FORMAT!'  LLBS.GE.60') 

1532  FORMAT! 15X, 'THETA ',7X,  '  SIGMA-SQUAR£',3X,' 
MERGING-DISTANCE:'/ 

1  3E20.8) 

1533  FORMAT!'  K,LLBS,TRACE!K ),THET A,DMRCE :'/2I5,3E20. 8) 
REWIND  NOFB 

REMIND  NUFC 

C  FOR  EACH  ITERATION:  ZERO  OUT  THE  VARIABLES 
DO  20  J=l,60 
CNO(J)=0. 

DO  10  L=l,2 
CKV!J,L )=0. 

10  CONTINUE 

20  CONTINUE 

LLBS=0 

C  DMIN:  SIGMA-SQUARE 

C  THETA:  SOME  THETA 

C  DMRGE:  MERGING  DISTANCE 

THETA=DMIN*DSTEP*FLOAT(IMRGE) 

DMRGE=S0RT!4./3.*!THETA-DMIN)) 

NR ITE! NOPF,l 532 )THETA,DMIN,DMRGE 
WRITE!?, 15G4HMRGE 

C  GO  THROUGH  SUBIMAGES  AND  LABEL  THEM  WITH  CLUSTERS 

DO  300  J=1,NQLH 
IF  (J.GT.l)  GOTO  35 

C  J=1  :  CASE  OF  FIRST  LINE  OF  SUBIMAGES 

DO  30  K=l,NOPH 
LABEL!2,K)=0 
DO  30  L  =1,2 
30  ABM!l,K,L)-0* 
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GOTO  45 

35  CONTINUE 

C  GET  PREVIOUS  LINE  OF  SUBIMAGES  IN  RGBM  ARRAY 

DO  40  K=l,NOPH 
DO  40  L=l/2 

40  ABM(1,K,L)=ABM(2,K,L) 

45  CONTINUE 

INDXB-J 

RE AD(NOFB"INDXB)( ABM(2,K, l),K=l,NOPH) 

INDXB= J+NOLH 

READ(NOFB'INDXB)(ABM(2,K,2),K=l,NOPH) 

C  INTIAL  LABEL  FOR  EACH  SUBIMAGE 

DO  50  K=l,NOPH 
LABEL(i,K)=LABEL(2,K) 

LABEL(2,K)=0 
50  CONTINUE 

INDXCsJ 

READ(NOFC'INDXC )( TRACE (K),K=1, NOP H) 

C  GO  THROUGH  IMAGES  ONE  BY  ONE 

DO  201  K=l,NOPH 

IF  (I0P(1).EQ.9)  WRITE(7,1533)K,LLBS,TRACE(K),THETA,DMRG£ 
IF  (LLBS.GE.60)  HRITE(7,1508) 

IF  (LLBS.GE.60)  GOTO  900 

C  CHECK  IF  THE  TRACE  OF  CURRENT  SUBIMAGE  >  THETA 

IF  (TRACE(K).GT. THETA)  GOTO  200 
C  SKIP 

IF  (J.GT.l)  GOTO  52 
C  J=l:  FIRST  LINE  OF  SUBIMAGES 

C  THE  FIRST  LINE  SECTION:  CONSIDERING  THE  NEIGHBOR 
Ml=2 
M2=2 
K1=K-1 
K2=K 
GOTO  54 

52  CONTINUE 

C  MOT  THE  FIRST  LINE;  SO  PREVIOUS  LINE  EXISTS 

Ml-1 
M2=2 
K1=K 
K2=K 

54  CONTINUE 

C  CHECK  IF  FIRST  SUBIMAGE  OR  NOT 

IF  (LL3S.EQ.0)  GOTO  90 
IF  (LABEL (Ml/K1)»EQ«0)  GOTO  55 

C  POTENTIAL  NEIGHBOR  NOT  LABELLED,  DIRECTLY  CHECK  CLUSTERS 

C  LA BEL( M1,K1 )  NEIGHBOR  HAS  BEEN  LABELLED 

C  AND  SPATIAL  CLUSTERING  SHOULD  BE  APPLIED 

DIFF=0. 

DO  62  L  =1,2 

62  DIFF=DIFF*(ABM(M1,K1,L)-ABM(M2,K2,L))**2 

DIFF=SQRT(DIFF) 

IF  (DIFF.GT.DMRGE)  GOTO  55 
C  WITHIN  MERGING  DISTANCE  ? 

LABEL ( M2, K2 )=L ABEL (Ml, Kl ) 

NL-LABEL(M1,K1 ) 
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LABEL! M2,K2 )=NL 
CNO(SL)=CNO(NL)+l. 

DO  64  L  =1/2 

64  CKV!NL,L)=!CKV!NL,L)*!CNO!NL)-l.)*ABM!M2,K2,L))/CNQ!NL) 
COTO  200 

55  CONTINUE 

DO  58  N=1,LLBS 
DIFF=0 • 

DO  56  L=l,2 

DIFF=DIFF*!CKV!N,L)-ABM!2,K2,L))**2 

56  CONTINUE 

FIRST! N)=SQRT!DIFF) 

58  CONTINUE 

CALL  DISMIN(FIRST,FMIN,NL,LLBS) 

IF  !FMIN.GT.DMRGE)  GOTO  90 
IF  !I0P!1).EQ.9)  WRITE!7, 1503)FMIN,DMRGE 
C  LABEL  CURRENT  SUBIMAGE  WITH  CLOSEST  CENTER 

LABEL! M2, K2)=NL 

C  UPDATE  NO.  OF  SUBIMAGES  OF  CURRENT  CLUSTER 

CNO!NL )=CNO!NL)+l. 

C  UPDATE  MEAN  VECTOR  OF  THIS  CLUSTER 

DO  60  L=i,2 

CKV!NL,L)=!CKV!NL,L)*!CNO!NL)-l.)+ABM!2,K2,L))/CNO!NL) 
60  CONTINUE 

IF  (I0P!1).EQ.9)  WRITE!?# 1506) 

GOTO  200 

C  NEW  CLUSTER  SECTION 

90  CONTINUE 

LLBS-LLBS+1 

C  UPDATE  OF  SUB  IMAGES  OF  THIS  CLUSTER 

CNO!LLBS)=CNO!LLBSm. 

C  UPDATE  NEW  CLUSTER  VECTOR 

DO  92  L=l,2 

92  CKV!LLBS,L)=ABM!M2,K,L) 

200  CONTINUE 

201  CONTINUE 

250  CONTINUE 

C  CHECK  CURRENT  LINE'S  LABELS 

IF  ! IOP (3 ).£Q. 1 )  WR ITE(7,1545 ) (LABEL (2,K),K=1,32) 

300  CONTINUE 

900  CONTINUE 

IF  !I0P!1).EQ.9)  WRITE! 7# 1502) (CNO!N),N=l,LLBS) 

IF  ! I0P!1).EQ. 9)  WRITE! 7# 1501 )IMRGE,LLBS 
WRITE! NOFF, 1501 )IMRGE,LLBS 
1545  FORMAT!'  LABEL', 3212) 

1501  FORMAT!/'  MERGE  ITERATION: ', 15, '  END  WITH  LLBS:  ',15/) 
RETURN 
END 
C 

C  SORTING  THE  KERNEL  VECTORS 

C 

SUBROUTINE  SORT!NCLRS) 

REAL  TEMP!2) 

COMMON  /BLOCKO/IQP!10),NOF1,NOF2,NOFB,NOFC,NOFF 
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COMMON  /BLOCK2/CKV(60,2),CNO(60),FKV(30,2),FNO(30) 

DO  30  I?2,NCLRS 
I2=NCLRS+1-I 
DO  20  J=l/I2 

IF  !CKV!J+1/1).GE.CKV!J,1))  GOTO  20 

DO  10  L  =1*2 

TEMP!L)*CKV!J*1/L) 

CKV(J+1,L)=CKV(J,L) 

CKV(J/L  )=TEMP(L) 

10  CONTINUE 

20  CONTINUE 

30  CONTINUE 

WRITE1NOFF/1533) 

1533  FORMAT!/*  SORTING  CANDIDATE  KERNEL  VECTORS*/) 

RETURN 

END 

C 

C  - 

C  TO  FIND  FINAL  KERNEL  VECTORS 

C  LIMIT  TO  10  ITERATIONS 

C 

SUBROUTINE  KERVEC(NC,KK,DIS) 

C  LIST  ARRAY  STORES  THE  TOTAL  DISTANCES  OF  ITERATIONS 

C  C  ARRAY  STORES  NUMBER  OF  PIXELS  FOR  EACH  CLUSTER 

C  D  ARRAY  STORES  TEMPORARY  DISTANCES  TO  CLUSTER  CENTERS 

C  FOR  CURRENT  PIXEL  BEING  PROCESSED 

REAL  DIST!10),FIRST!60) 

COMMON  /BLOCKO/IOP(10)#NOF1,NQF2#NQFB,NQFC,NOFF 
COMMON  /BLOCK1/NOL,NOP,NOL2,NOP2,NOLH,NOPH 
COMMON  /BLOCK2/CKV(60,2),CND(60),FKV(30,2),FNO(30) 
COMMON  /BLOCK3/IA(256),IB(256),A(256,2),B<256,2> 
COMMON  /BLOCK 4/ ABC256, 2)/ ABM( 2,128,2), TRACEC128) 
COMMON  /BLOCK5/LABEL(2,128),ICOLOR(256),IBtf(256) 

1020  FORMAT! *  IN  KERVEC#  KKs*,I5) 

C  FINAL  KERNELVECTORS  SAVED  IN  FKV  ARRAY 
C  FNO  STORE  NO.  OF  PIXELS  IN  EACH  CLUSTER 

C  K-MEANS  ALGORITHM 

C  10-JUN-82  CORRECT  IMPLEMENTATION 

C  REFERENCE:  TOU  AND  GONZALEZ/11  PATTERN  RECOGNITION 
C  PRINCIPLES"/  PP.94-97/ 11974). 

DO  12  N=1/NC 
FNO!N)=CNO!N) 

CNO(N)=0. 

DO  10  L-1,2 
FKV1N, L)=CKV! N/L) 

CKV!N,L  )=0. 

10  CONTINUE 

12  CONTINUE 

C  NO.  OF  ITERATIONS  LIMIT  TO  10 

DO  500  KK=1/10 
ERITE17/  1020 )KK 
NR ITE! N OFF/ 1020 )KK 
DO  15  N=1,NC 
15  FNO!N)=0. 

REWIND  NOF1 
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REWIND  NOF2 
RE HIND  MOFB 

C  FOR  EACH  ITERATION: 

C  REWIND  THE  FEATURES  FILES:  A  ,  B  COMPONENTS 
C  REHIND  THE  TEMPORARY  CLASSIFIED  RESULT  FILE/  NOFB 
C  CLASSIFYING  STANDARDS  IN  FKV  ARRAY 

C  AT  THE  SAME  TIME/  COLLECTING  THE  NEH  CENTERS  IN  CKV  ARRAY 
C  I.E./  UPDATING  THE  KERNEL  VECTORS  BY  CURRENT  CLUSTERING 
C  ICOLOR  ARRAY  STORES  CLASSIFIED  RESULT  OF  CURRENT  LINE 
DO  200  1=1/ NOL 
INDX1=I 

READ!N0F1'INDX1)!IA!K),K=1,N0P) 

1NDX2=I 

READ!NOF2*INDX2)!IB!K),K=l,NOP) 

DO  20  J=l,NOP 
AB(J/1 ) =FLOAT!IA! J) ) 

AB(J/2 )=FLOAT(IB(J)) 

20  CONTINUE 

C  GO  THROUGH  PIXELS  TO  LABEL  THEM  HITH  CLUSTERS 
DO  100  J=l,NOP 
DO  40  N=1,NC 
SUM=0. 

DO  30  L=l/2 

30  SUM=SUM+!AB!J,L)-FKV!N,L))**2 

40  FIRST! N )=SQRT(SUM) 

CALL  D ISMIN(FIRST,FMIN/NN/NC> 

ICOLOR  !J)=NN 
CNO!NN)=CNO!NN)+l. 

C  CURRENT  PIXEL  HAS  FOUND  CLOSER  TO  NN-TH  CLUSTER 

C  THE  FEATURES  SHOULD  BE  INCLUDED  TO  UPDATE  THE  NN-TH 
C  KERNEL  VECTOR 

DO  70  L  =1/2 

70  CKV!NN,L)=!CKV!NN,L)*!CNO!NN)-l.)*AB!J,L))/CNO!NN) 

100  CONTINUE 

INOXB=I 

IF  (I3P<2).EQ.l)  WRITE(7/1505)( ICOLOR (K)/K =1,64) 

WRITE! NOFB* I NDXB)( 1 COLOR! K),K=1, NOP) 

200  CONTINUE 

C  CURRENT  IN  CKV)  STORE  IN  FKV  TO  BE  USED  TO  CLASSIFY 
C  IN  ROUTINE  D1STAN,  AND  HILL  GIVE  TOTAL  DISTANCE 

DO  260  N=1,NC 
FNO!N)=CNO!N) 

CNO!N)=0. 

DO  250  L=l,2 
FK V ! N, L )=CKV ! N/L) 

CKV!N/L)=0. 

250  CONTINUE 

260  CONTINUE 

WRITE! NOFF,1503)!!FKV!N/L)/L=1/2),N=1,NC) 

1503  FORMAT! /10X, 'CHECK  FKV:  */ 

130!!5X,2E20. 8)/)) 

WRITE! NOFF/1506)!FNO!N),N=1, NC) 

1506  FORMAT!*  OF  PIXELS  IN  EACH  CLUSTER:*/ 

1  30!/Fl2.1)) 

C  ******** 
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CALL  DISTAN!NC,DIS) 

C  ******** 

DIST(KK)=DIS 

WRITE!7,1504)DIST!KK),KK 
WRITE! NOFF,1504)DIST!KK),KK 

1504  FORMAT!/'  IN  KERVEC:  DIST!KK)=',E20. 8, '  KK=',I2/) 
IF!KK.EQ.1)G0  TO  500 

RATIO* ! DI ST!KK )~DIST!KK-i ))/D 1ST! KK-1) 

WRITE! NOFF,1005)RATIO 
IF (ABS! RATIO).LT. 0.001)  GOTO  900 
500  CONTINUE 
900  CONTINUE 

1005  FORMAT!'  RATIO  IN  KERfECS *, E20. 8) 

1505  FORMAT!'  COLORS ',6411) 

RETURN 

END 

C 

C  USE  CURRENT  KERNEL  VECTORS  TO  CALCULATE  THE  TOTAL 
C  DISTANCE  OF  THE  IMAGE 

C  NCS  NUMBER  OF  CLUSTERS 

C  DISTOTS  !  RESULT  )  TOTAL  DISTANCE 
C 

SUBROUTINE  DISTAN!NC,DISTOT) 

COMMON  /BLOCKO /  IOPUO), N0F1 , N0F2, NOF B, NOFC, NOFF 
COMMON  /BLQCK1/NOL,NOP,NOL2,NOP2,NOLH,NOPH 
COMMON  /BLOCK2/CKV!60,2 ),CNO!60),FKV !30,2),FNO!30) 
COMMON  /BL0CK3/IA!256),IB!256),A!256,2),B!256,2) 
COMMON  /BL0CK4/AB!256,2),ABM!2,128,2),TRACE!126) 
COMMON  /BL0CK5/LABEL!2,128),IC0L0R!256),IBW!256) 
DISTOTsO. 

REWIND  N0F1 
REWIND  N0F2 
REWIND  NOFB 
DO  200  Isl,NOL 
INDXB=I 

READ!NOFB'INDXB)!ICOLOR!K),K=1,NOP) 

INDX1=I 

READ!N0F1'1NDX1)(IA(K),K=1,N0P) 

INDX2SI 

RE AD!NQF2 'INDX2 )! IB!K),K*1, NOP) 

C  STORE  FEATURE  VECTOR  IN  WORKING  VARIABLE  X 

DO  10  K=1,N0P 
AB!K,1 )=FLOAT(IA(K)) 

AB!K,2 )=FLOAT(IB(K ) ) 

10  CONTINUE 

DO  100  J=1,N0P 
NCLSRsICOLOR! J) 

SUM=0. 

DO  30  L  =1,2 

SUM=SUM+( AB (J,L)-FKV(NCLSR,L) ) **2 
30  CONTINUE 

D I STOTsDI STOT+SQRT ( SUM ) 

100  CONTINUE 

200  CONTINUE 
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NR ITE(NOFF/ 1501)DI STOT 

1501  FORMAT!/'  IN  DISTAN:  OISTOT  =  ',E20.8/) 

RETURN 
END 

C 
C 

C  USE  THE  FINAL  KERNEL  VECTORS  TO  CLASSIFY  THE  IMAGE 

C  INTO  CLUSTERS  (  SUBREGIONS  ) 

C 

SUBROUTINE  CLASS(NC) 

REAL  FIRSTC60) 

COMMON  /BLOCK 0/  IOP  ( 1 0  )  ,  NOF 1,  N OF2,  NOF  B,  NOFC/  NOFF 
COMMON  /BLOCK1/NOL/NOP/NQL2/  NOP2/NOLH/NOPH 
COMMON  /BLOCK2/CKV(60,2),CNO(60),FKV(30,2),FNO(30) 
COMMON  /BLOCK3/IA!256)/IB(256)/A(256/2)/B(256/2) 

COMMON  /BLOCK 4/ AB  (256/2  )/ ABM  (2#  128/2  )/ TRACE  (128) 

COMMON  /BL0CK5/LABEL!2/128)/IC0L0R!256),IBW(256) 

WRITE! NOFF# 1505 ) 

1505  FORMAT!/'  IN  CLASS  S'/) 

WRITE(7/1501)((FKV(N/L),L=1,2),N=1/NC) 

WRITE! NOPF, 1501 )((FKV(N/L),L=1,2),N-1,NC) 

1501  FORMAT! 10X/'  FINAL  KERNEL  VECTERS  : '/(30(5X/2E20 .8/ ) ) ) 
1509  FORMAT! 1X/64I1) 

C  USE  FINAL  KERNAL  VECTORS  TO  CLASSIFY  THE  PICTURE 
C  NOFC:  INTEGER  OUTPUT;  NOFD:  REAL  OUTPUT 

REWIND  NOF1 
REWIND  NOF2 
REWIND  NOFB 
REWIND  NOFC 

C  RANGE  OF  FINAL  RESULT  0..180 
F  ACT=180. /FLOAT(NC) 

DO  200  I=l/NOL 

INDX1-I 

INDX2=I 

READ(NOFl'INDXl)(IA!K),K=l,NOP) 

READ(NOF2 'INDX2 )(IB(K)/K=l/NOP) 

DO  10  J=l/MOP 
AB(J/1 ) -FLOAT! I A! J) ) 

AB(J/2 )=FLOAT (IB(J>) 

10  CONTINUE 

DO  100  J=l/NOP 
DO  40  N=1,NC 
SUM=0. 

DO  30  L  =1/2 

SUM-SUM*! AB(J/L)-FKV(N/L) )**2 
30  CONTINUE 

FIRST! N  )  =  SORT (SUM) 

40  CONTINUE 

C  DISTANCES  TO  FINAL  KERNEL  VECTOS  OF  NC  CLUSTERS 

C  FROM  CURRENT  PIXEL  ARE  STORED  IN  DIS  ARRAY/ 

C  CALLING  SUBROUTINE  DISMIN  TO  FIND  TO  WHICH  CLUSTER 
C  THE  CURRENT  PIXEL  IS  CLOSER  (  MINIMUM  DISTANCE  ). 

C  THE  RESULT  IS  KMIN-TH  CLUSTER 

CALL  DISMIN(FIRST,SMIN,NMIN/NC) 


o  « 
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BLACK  AND  WHITE  DISPLAY  PURPOSE:  NEEDS  TO  MULTIPLY  A  PACTOR 
TO  BE  IN  REASONABLE  GRAY  LEVEL  RANGE 
IBW(J)sINT<PLOAT(NMIN)*FACT) 

C  COLOR  DISPLAY  PURPOSE:  AN  INTEGER  IN  RANGE  1  TO  7 
ICOLOR<J)xNMIN 
100  CONTINUE 
1NDXB=I 

WRITE(MQFB'INDXB)(IBW(K)*K-l*NOP) 

INOXCsI 

NRITE(NOFC'INDXC)(ICQLOR(K)*K*l*NOP) 

IP  (I.LE.64)  WRITE(7,1509)(ICOLOR(K)*K=1*64) 

IF  (I.LE.64)  NRITE(NOPF*1509)(ICQLQR (K)*K=1*64) 

200  CONTINUE 
RETURN 
END 
C 

C 

SUBROUTINE  DISMIN(DARRY*DATMIN*NOMIN*NCLSTR) 

C  PASS  DARRY  ARRAY  WITH  NCLSTR  ELEMENTS  MEANINGFUL 
C  SEARCH  FOR  THE  MINIMUM  ELEMENT*  NOMIN-TH  ELEMENT* 

C  WITH  VALUE  DATMIN;  PASS  BACK  DATMIN  AND  NOMIN  BACK 
C  CALLING  ROUTINE 

REAL  DARRYC60) 

DATMIN=DARRY(1) 

NOMIN-1 

C  ASSUME  THE  FIRST  ELEMENT  IS  THE  MINIMUM 

C  THEN  GO  THROUGH  THE  REST  OP  THE  ARRAY  TO  FIND  AMY  SMALLER 

IF  (NCLSTR. EG.  1)  GOTO  900 
DO  100  1-2* NCLSTR 
IF  (DARRY(I).GE. DATMIN)  GOTO  100 
DATMIN-DARRY(I) 

NOMIN- I 

100  CONTINUE 
900  CONTINUE 
RETURN 


APPENDIX  SPACLR.FOR 
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9C12 

C  FILE  NAME:  SPACLR.FOR 

C  AUTO  SPATIAL  CLUSTERING  FOR  COLOR  IMAGE  DATA 

C 

C  REFERENCE  FUN ADA/"  SPATIAL  CLUSTERING  PROCEDURESFOR  REGION 
C  ANALYSIS"/ PATTERN  RECOGNITION, 1 2,395-403, ( 1980). 

C 

C 

C  INPUTS:  NOF1,  FTN1.DAT,  FEATURE  1 

C  NOF2,  FTN2.DAT,  FEATURE  2 

C  NOF3,  FTN3.DAT,  FEATURE  3 

C  OUTPUTS :NOFC,  FTN12.DAT,  CLUSTERING  RESULT  FOR  COLOR 
DISPLAY 

C  NOFD,  FTN11.DAT/  CLUSTERING  RESULT  FOR  BLACK/HHITE 

C  DISPLAY 

C  NOFF/  FTN15.DAT,  SOME  PARAMETERS  DURONG  PROCESSING 

C  SIZE  OF  IMAGE  IS  RESTICTED  TO  256  BY  256  IN  FEATURE  SPACE 


LOGICAL*l  DDMMYC9) 

REAL  PREY (60,3 )/PRENO( 60) 

COMMON  /BLQCKO/IOP(10)/NOFl,NOF2,NOF3,NOFC,NOFD,NOFF 

COMMON  /BLOCK 1/NOL/ NOP/ NOL2/NOP2/NOLH, NOP H 

COMMON  /BLOCK2/CXV(60/3)/CNO(60),FKY(30,3),FNO(30) 

COMMON  /BLOCKS/ IR (256)# 16(256), IB (256), RGB (256, 3) 

COMMON  /BL0CK4/RGBM(2,128,3)/TRACE(128) 

COMMON  /BLOCK5/LABEL(2,128),ICOLOR(256),IBH(256) 

COMMON  /BLOCK6/R(256,2)/G(256,2)/B(256,2) 

DATA  NOF1/NOF2/ NOF3, NOFC/ NOFD, NOFF /l/ 2/ 3, 12, 11, 15/ 

1701  FQRMAT(*  RUNNING  PROGRAM  SPACLR.FOR  AUTO  SPATIAL 
CLUSTERING*/ 

1  *  INPUTS  ,  OUTPUTS:  NOF1, NOF2, NOF3, NOFC, NOFD, NOFF '/ 

6  *  FOR  CHECK:', 615/ 

2  '  NOFC  ,  NOFD  STORE  TEMPORARY  DATA  DURING 
PROCESSING'/ 

3  '  OUTPUTS:  NOFC(REHRITTEN)  7  OR  LESS  COLORS'/ 

4  '  NOFD (REWRITTEN)  BLACK  AND  WHITE 

DISPLAY'/ 

5  '  NOFF,  INFORMATIONS  DURING  PROCESSING') 

1515  FORMAT!'  MERGE  ITER ATION', 15 ) 

1520  FORMAT! '  THE  ',I5,'-TH  ITERATION  REACHES  MAXIMUM  NO. 
CLUSTERS'/ 

1  '  NO.  OF  CLUSTERS:', 15) 

1511  FORMAT!'  KERNEL  CANDIDATE  VECTORS  FOR', 14,'  STARTING  ', 

1  'CLUSTER  CENTERS'/) 

1512  FORMAT!*  ENTER  IMAGE  DATA  FILE  SIZE  (  FEATURE  SPACE  )'/ 

1  '  NOL,  NOP:  F0RMAT(2I5)') 

1513  FORMAT(2I5) 

1514  FORMAT!'  CHECK:NOL,NOP,NOL2,NOP2,NOLH,NOPH*/6I5) 

1516  FORMAT!'  ENTER  OPTIONS:  ( IOP(K),K=l, 10) '/ 

1  '  IOP(l):  CONTROLS  PRINTER,  1:  MEANS  ,  TRACE 

MATRIX*/ 
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2  *  IOPC2):  K-MEANS  ALGORITHM  DETAILS  ON  SCREEN*/ 

3  *  IQPC3):  MERGE  DETAILS#  LABEL (2#K)  ARRAY*/ 

4  *  I0PC4):  NSTEP#  NO.  OF  STEPS*/ 

5  '  I0PC5):  1,  SKIP  K-MEAN  ITERATION*/) 

1517  FORMAT (10 16) 

1518  FORMAT ( *  TODAY  IS  *#9A1) 

REMIND  NOFF 

MR ITE( 7*1501) 

MR ITE( NOFF# 1501) 

CALL  OATE(DDMMY) 

MRITEC  NOFF# 1518 )(DDMM¥(K)#K=1# 9) 

MRITEC  7/1510)  (DDMMYC  K)#K-1#  9) 

MRITEC  M OFF# 1701 )NOFl#NOF2#NOF3#NQFC#NOF D# NOFF 
MR 1TE( 7#1701 )NOFl#NOF2#NOF3#  NOFC#NOFD#NOFF 
MRITE(7#1512) 

MRITEC NOFF#1512) 

RE AD (5# 1513) MOL# NOP 
MRITE(7,1516) 

MRITEC  NOFF# 1516) 

READC5# 1517)(IOP(K)#K=1#10) 

MRITEC NOFF# 1517)(IOP(K)#K=l# 10) 

MRITEC  7#1517)C10PCK)#Ksl#10 ) 

NOL2-NQL+NOL 

NOP2=NOP*NOP 

NOLH-NOL/2 

NOPH=NOP/2 

MRITEC  NOFF#1514 )NOL#NOP#NOL2# NOP2,NOLH, NOPH 
HR  IT EC  7 # 1 514 ) NOL# HOP #NOL2# NOP 2 #NOLH# NOP H 
DEFINE  FILE  N0F1CN0L#N0P#U#INDX1) 

DEFINE  FILE  NOF2(NOL#NOP#U,INDX2) 

DEFINE  FILE  NOF3(NOL#NOP#D# INDX3) 

DEFINE  FILE  NOFC(NOL#NOP#U# INDXC) 

DEFINE  FILE  NOFD(NOL#NOP#U#INDXD) 

C  N0FF=FTN15.DAT  UNFORMATTED 
REMIND  NOF1 
REMIND  NOF2 
REMIND  NOF3 
REMIND  NOFC 
REMIND  NOFD 

1501  FORMAT! '  THIS  IS  THE  LOG  FILE  OF  EXECUTING  SPACLR.FOR*) 

C  CALCULATE  MEANS  OF  FEATURE  VECTORS  OF  2  BY  2  SUB  IMAGE 
C  STORE  IN  NOFC  AND  HALF  NOFD  128  BY  128  SUBIMAGES 

C  IN  NOFC:  FIRST  HALF  IS  R  MEANS#  SECOND  HALF  IS  G  MEANS 

C  IN  NOFD:  FIRST  HALF  IS  B  MEANS#  SECOND  HALF  IS  TRACES 

CALL  DISPER 


C  FIND  MAX  #  MIN  OF  TRACE  MATRIX 
CALL  MAXMIN(DMAX#DMIN) 


C  MERGING  SECTION 
NSTEP=IQP(4) 

STEP=( DMAX-DM IN )/PLOAT( NSTEP) 
MRITEC 7 #1522 ) DM AX# DM IN# STEP 
MRITEC NOFF#1522)DMAX#DMIN#STEP 
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1522  FORMAT!  *  DMAX**,E20.8,*  DMIN=*,E20.8,*  STEP**,E20.8) 
IPREfsO 

C  ITERATIONS  TO  FIND  MAXIMUM  NO.  OF  CLUSTERS 
DO  300  I»1,NSTEP 
IM=I 

URITE(NOFF,1515)IM 

CALL  MERCE!IM,STEP,DMIN,NCLSR) 

C  ACCEPTED  NO.  OF  CLUSTERS:  7  OR  LESS 

IF  (IPREV.LE.7. AND. IPREV.GT.l )  GOTO  333 
IPREV-MCLSR 

C  SAVE  CURRENT  NUMBER  OF  CLUSTERS  AND  KERNEL  VECTORS 
DO  200  Jsl,IPREV 
PRENO( J )=CNO( J) 

PREV( J* 1)=CKV( J,l) 

PREV(J,2)=CKV(J,2) 

PREV(J#3)sCKV(J/3) 

200  CONTINUE 
300  CONTINUE 
333  NI=IM-1 

HR ITE( HOFF/1520 )NI, IPREV 
DO  350  J=1,IPR£V 
CN0( J) =PRENO( J) 

CKV(J# 1 )~PREV(Jsl) 

CKV(J,2)=PREV(J,2) 

CKV( J/3 )sPREV(J/3) 

350  CONTINUE 

HR ITE( MOFF/1561 ) 

1561  FORMAT! •  MERGE  ENDEO  WITH  MAXIMUM  NO.  CLUSTERS') 
WRITE!  HOFF#  1562)!!  CKV!N/L)#L-1#3)^N=1#  IPREV) 

1562  FORMAT!*  BEFORE  SORTING'/ 

1  *  KERNEL  CANDIDATE  VECTORS*/ 

2  30(!5X,3E20.8)/)) 

C  SORT  THE  CANDIDATE  VECTORS 

NC=IPREV 


SORT  THE  CANDIDATE  KERNEL  VECTORS 
CALL  SORT!NC) 

NRITE(NOFF,1563>!!CKV<N,L),L=l,3),N=l,NC) 

FORMAT!*  SORTED  KERNEL  CANDIDATE  VECTORS'/ 

1  30((5X,3E20.8)/)) 

IF  (NC.GT.7)  NC=7 

FOR  THE  PURPOSE  OF  AED-512  PSEUDO  COLOR  DISPLAY 
URITE!N0FF,1511)NC 

IF  !IOP!5).EQ.l)  SKIP  THE  K-MEAN  ITERATIONS 
DIRECTLY  USE  MERGING  RESULT  CADIOATE  KERNEL  VECTORS 
TO  CLASSIFY  THE  IMAGE 
HR1TE!7,1568)IQP!5) 

FORMAT!*  IOP(5)s*,I5) 

IF  !I0P!5).NE.l)  GOTO  700 
HR ITE!7/1570) 

FORMAT!*  SKIP  K-MEAN  ITERATION*) 

DO  650  N*1,NC 
DO  640  L-1,3 
FKV(N/L)-CKV(N/L) 

CONTINUE 
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650 

700 

1580 

C 

C 

C 


1500 

800 

C 

C 

C 

C 


1523 


C 

C 

C 

c 

c 

c 

c 


c 


c 

10 

20 


CONTINUE 
GOTO  800 
CONTINUE 
MRITE(7,1580) 

FORMAT { *  CALLING  KERVECt  K-MEAN  ITERATION') 


ITERATIONS  TO  FIND  MORE  ACCURATE  KERNEL  VECTORS 
CALLED  FINAL  KERNEL  VECTORS 
CALL  KERVEC(NC,KK,DD) 

URITE(NOFF,1500)KK,DD 

FORMAT (11/ 'CLUSTERING  REPEATS', IX, 13, IX, 'TIMES'/IX, 
l'THE  FINAL  UITHIN-CLASS  DISTANCE  IS*,1X,E20.8/) 
CONTINUE 


CLASSIFICATION  SECTION 

OUTPUTS:  NOFC,  COLOR  DISPLAY  RESULT 

NOFD,  BLACK/WHITE  DISPLAY  RESULT 
CALL  CLASS(NC) 

HR ITE( HOFF# 1523) 

FORMAT ( 10 X, ' 1 1 1  COMPLETE  EXECUTION  OF  PROGRAM  SPACLR  III') 

CALL  EXIT 

END 


-  SUBPROGRAMS  - 


SUBROUTINE  TO  CALCULATE  TRACE  MATRICES  OF  FEATURE  MATRICES 
STORED  IN  NOFD 
SUBROUTINE  DISPER 

COMMON  /aLOCKO/IOPUO),NOFl,NOF2,NOF3,NOFC,NOFD,NOFF 
COMMON  /BLOCKl/NOL,NOP,llOL2,NOP2,NOLH,NOPH 
COMMON  /BLOCK3/IRC256), IGC256), IB(256), RGB(256,3) 

COMMON  /BLOCK4/RGBM(2,128,3 ),TRACE(128) 

COMMON  /BLOCK6/R(256,2),G(256,2),B(256,2) 

REWIND  NOF1 
REMIND  NOF2 
REMIND  NOF3 
REMIND  NOFC 
REHIND  NOFD 

PROCESS  THROUGH  ROUS  OF  DATA  MATRIX 
DO  100  I=l,NOLH 
12=1*1-1 
INDX1=I2 
INDX2=I2 
INDX3= 12 
DO  40  JJ=1,2 

READ  2  LINES  OF  EACH  FILE 
READ(N0F1 'INDX1 )( IR (K), K=l, NOP ) 

DO  10  J=l,NOP 

R( J, JJ )=FLOAT(IR(J)) 

READ(NOF2'INDX2)(IG(K),K=l,NOP) 

DO  20  J=l,NOP 

G( J, JJ ) =FLOAT(IG( J) ) 

READ(NOF3'INDX3  )(IB(K),K=l,NOP) 
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DO  30  J=l,NOP 
30  B(J,JJ)=FLOAT(IB(J)) 

40  CONTINUE 

C  CALCULATION  THROUGH  EACH  SUB  IMAGE 

DO  80  Ksl#NOPH 
K1=IUK-1 
K2=KU1 
S1=0. 

S2=0. 

S3=0. 

DO  62  M=K1,K2 
DO  60  L  =1/2 
S1=S1+R(M,L) 

S2=S2*G(M/L) 

S3=S3+B!M,L) 

60  CONTINUE 
62  CONTINUE 

RGBM(2,K,1)=S1*0.25 

RGBM(2,K<,2)=S2*0.25 

RGBM(2,K,3)=S3*0.25 

S1=0. 

S2=0. 

S3=0  . 

DO  72  M-K1,K2 
DO  70  L=l#2 

S1=S1*(R(M,L>-RGBM(2,K,1))**2 
S2=S2«-(G(M,L )-RGBM(2#K*  2))**2 
S3-S3+! B( M,L)-RGBM(2,K,3))**2 
70  CONTINUE 

72  CONTINUE 

TRACECK)=(S1*S2*S3)*0.25 
80  CONTINUE 

INDXC=I 

IF  (IOP(l).EQ.l)  WRITE(6,1001)(RGBM(2,K,1),K=1,32) 

NRITE(NOFC'INDXC)(RGBM(2,K,l),K=l,NOPH) 

INDXC=I*NOLH 

IF  (IOP(l).EQ.l)  NRITE(6,1002)(RGBM(2,K,2),K=1,32) 
NRITE(NOFC'INDXC)(RGBM(2#K#2)#K=l,NOPH) 

INCXDsI 

IF  (IOP(l).EQ.l)  WRITE!  6#  1003) (RGBM( 2/K/3)sK=ls32) 
WRITE!  NOFD'INDXD)  (RGBM(  2#K#  3  )/K=l*  NOPH) 
INDXDsI+NOLH 

IF  (lOP(l).EQ.l)  WRITE(6,1004)(TRACE!K>,K=1,32) 
WRITE!  NOFD'INDXD)  (TRACE  (X)/ K=l*NOPH) 

100  CONTINUE 

1001  FORMAT! *  RM',32F4.0) 

1002  FORMAT!'  GM',32F4.0) 

1003  FORMAT!*  BM',32F4.0) 

1004  FORMAT!'  TR',32F4.0> 

RETURN 

END 

C  - 

C  READ  TRACE  MATRIX  TO  FIND  MAX  ,  MIN 
C 

SUBROUTINE  MAXMIN(DMAX, DMIN ) 
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I 


1 


COMMON  /BLOCKO/ IOP ! 10 ) , NOF1, NOF2, NOF3, N OFC, NOPO, NOFF 
COMMON  /BLOCK1/NOL, NOP, NOL2, NOP2, NOLH, NOPH 
COMMON  /BLOCK4/ RGBM( 2, 1 28,3 ), TRAC E( 128) 

REMIND  NOFD 
INDXD=l+NOLH 

READINOFD'INDXD )! TRACE (K),K= 1, NOPH) 

DMAX=TR ACE(l) 

DMIN=TRACE!1) 

DO  10  J=2,NOPH 

IF  (TRACE(J).LT.DMIN)  DMIN=TRACE(J) 

IF  (TRACE(J).GT.DMAX)  DMAX=TR ACEC J) 

10  CONTINUE 

DO  100  1-2, NOLH 
1NDXD=I+N OLH 

READ!NOFD *INDXD)!TR ACE! K),K=1, NOPH) 

DO  30  J=l,NOPH 

IF  !TRACE!J).LT.DMIN)  DMIN=TRACE( J) 

IF  (TRACE(J).GT.DMAX)  DMAX=TRACE( J) 

30  CONTINUE 
100  CONTINUE 

MR IT£( 7, 1001) DM AX, DM IN 

1001  FORMAT!/*  DMAX=*,F12.4,*  DMIN=*,F12. 4) 

RETURN 

END 

C 

C  - 

C  MERGE  AND  DECIDE  KERNEL  CANDIDATE  VECTORS 

C 

SUBROUTINE  MERGEC IMRGE,DSTEP, DMIN,LLBS) 

REAL  FIRSTC60) 

COMMON  /BLOCKO/ IOP !10),NOF1,NQF2,NOF3,NOFC, NOFD, NOFF 
COMMON  /BL0CK1/N0L,NQP,N0L2, NOP2, NOLH, NOPH 
COMMON  /BLOCK2/CKV (60,3 ),CN 0(60 ),FKV (30,3),FNO(30) 
COMMON  /BLQCK4/RGBM(2,128,3 ), TRAC E( 128 ) 

COMMON  /BLOCK5/LABEL(2, 128), ICOLOR(256),IBM(256) 

1502  FORMAT(30(/'  CNO! N)  :  *,F12. 1 )  ) 

1503  FORMAT!*  *FMIN,  DMRGE: *,2E20 • 8 ) 

1504  FORMAT!/*  IMRGES*,I5) 

1505  FORMAT!*  FIRST  SUBIMAGE:') 

1506  FORMAT!*  JOIN  CLUSTER  NL:*) 

1507  FORMAT!*  NEM  CLUSTER:*) 

1508  FORMAT!*  LLBS.GE.60*> 

1532  FORMAT! 15X, 'THETA*, 7X, *  SIGMA-SQUARE*, 3X,* 
MERGING-DISTANCE:*/ 

1  3E20.8) 

1533  FORMAT!*  K,LLBS,TRACE!K ),THET A,DMRGE :*/2I5,3E20. 8) 

REMIND  NOFC 

REMIND  NOFD 

C  FOR  EACH  ITERATION:  ZERO  OUT  THE  VARIABLES 

DO  20  J=l,60 
CNO(J)=0. 

DO  10  L  =1,3 
CKV!J,L)-0. 

10  CONTINUE 

20  CONTINUE 


I 
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LLBS-0 

C  DMIN:  SIGMA-SQUARE 

c  theta:  some  theta 

c  dmrge:  merging  distance 

THETA=DMIN+DSTEP*FLOAT! IMRGE) 

DMRGE=SQRT!4./3.*! THET A-DMIN) ) 

WRITE!  N OFF/ 1532 )TH£TA/ DMIN* DMRGE 
WRITE! 7/1504) IMRGE 

C  GO  THROUGH  SUBIMAGES  AND  LABEL  THEM  WITH  CLUSTERS 

DO  300  J=l/NOLH 
IF  (J.GT.l)  GOTO  35 

C  J=1  I  CASE  OF  FIRST  LINE  OF  SUBIMAGES 

DO  30  K=l/NOPH 
LABEL! 2/K)=0 
DO  30  L=l/3 
30  RGBM(1,K,L)=0. 

GOTO  45 

35  CONTINUE 

C  GET  PREVIOUS  LINE  OF  SUBIMAGES  IN  RGBM  ARRAY 

DO  40  K-l/NOPH 
DO  40  L  =1/3 

40  RGBM(1,K,L)=RGBM(2,K,L) 

45  CONTINUE 

INDXC=J 

RE AD (NOFC 'INDXC )(RGBM( 2/K/l )/K=l#NQP H) 

INDXC= J+NOLH 

READ(NQFC"INDXC)(RGBM(2/K/2 )/K=l#NOPH) 

INDXD-J 

READ!NOFD'INDXD)!RGBM!2/K/3)/K=l/NOPH) 

C  INTI AL  LABEL  FOR  EACH  SUBIMAGE 

DO  50  K=l/NOPH 
LABEL!1/K)=LABEL!2/K) 

LABEL! 2/K)=0 
50  CONTINUE 

1NDXD=J*N0LH 

READ!NOFD'INDXD)!TRACE!K)/K=l/NOPH) 

C  GO  THROUGH  IMAGES  ONE  BY  ONE 

DO  201  K=l#NOPH 

IF  !  IIP !1 ) . EQ. 9  )  WRITE!7/1533 )K/LLB S/TRACE !k)/ THETA/ DMRGE 
IF  !LLBS.GE.60)  WRITE! 7/1508) 

IF  1LLBS.GE.60)  GOTO  900 

C  CHECK  IF  THE  TRACE  OF  CURRENT  SUBIMAGE  >  THETA 
IF  !TRACE!K).GT. THETA)  GOTO  200 

C  SKIP 

IF  ! J. GTa 1 )  GOTO  52 
C  J=l:  FIRST  LINE  OF  SUBIMAGES 

C  THE  FIRST  LINE  SECTION:  CONSIDERING  THE  NEIGHBOR 

Ml=2 
M2=2 
K1=K-1 
K2=K 
GOTO  54 

52  CONTINUE 

C  NOT  THE  FIRST  LINE;  SO  PREVIOUS  LINE  EXISTS 

Ml  =  l 
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M2=2 

K1=K 

K2=K 

54  CONTINUE 

C  CHECK  IF  FIRST  SUBIMAGE  OR  NOT 

IF  (LLBS.EQ.O)  GOTO  90 
\  IF  !LABEL!Ml*Kl).EQ.O)  GOTO  55 

q  POTENTIAL  NEIGHBOR  NOT  LABELLED*  DIRECTLY  CHECK  CLUSTERS 

C'<  LABEL(M1,K1)  NEIGHBOR  HAS  BEEN  LABELLED 

C  AND  SPATIAL  CLUSTERING  SHOULD  BE  APPLIED 

DIFF=0. 

DO  62  L»l*3 

62  DIFF=DIFF+(RGBM(M1*K1,L)-RGBM(M2*K2,L))**2 

DIFF=SGRT(DIFF) 

IF  CCIFF.GT.DMRGE)  GOTO  55 
C  WITHIN  MERGING  DISTANCE  ? 

♦  LABEL! M2* K2 )=LA8EL! M1*K1) 

NL-LABEL(M1*K1) 

LABEL!  M2* K2  )=NL 
CNO!NL)=CNO!NL)+1. 

DO  64  L=l*3 

64  CKV!NL*L)=!CKV!NL*L )*! CNO! NL)-1. J+RGBM! M2,K2*L))/CNO!NL) 

GOTO  200 

55  CONTINUE 

DO  58  N=1*LLBS 
DIFF=0. 

DO  56  L=l*3 

DIFF=DIFF+!CKV!N*L)-RGBM<2*K2*L))**2 

56  CONTINUE 

FI RST! N  )=SQRT!DIFF) 

58  CONTINUE 

CALL  DISMIN!FIRST*FMIN*NL*LLBS) 

IF  ! FMI N« GT. DMRGE)  GOTO  90 
IF  !I0P!1).EQ.9)  MRI TE ! 7*1503 )FM IN* DMRGE 
C  LABEL  CURRENT  SUBIMAGE  WITH  CLOSEST  CENTER 

LABEL!  M2*  K2)=NL 

C  UPDATE  NO.  OF  SUBIMAGES  OF  CURRENT  CLUSTER 

CNQ!Nu)=€NQ!NL)4-l. 

C  UPDATE  MEAN  VECTOR  OF  THIS  CLUSTER 

DO  60  L=1 *3 

CKV!NL*L)=!CKV!NL,L)*!CNQ!NL)-1.)*RGBK!2,K2*L))/CNQ!NL) 
60  CONTINUE 

IF  !I0P!1).EQ.9)  WRITE! 7*1506) 

GOTO  200 

C  NEW  CLUSTER  SECTION 

90  CONTINUE 

LLBSsLLBS+1 

C  UPDATE  OF  SUBIMAGES  OF  THIS  CLUSTER 

CNC!LLBS)=CNO!LLRS )+l. 

C  UPDATE  NEW  CLUSTER  VECTOR 

DO  92  L=l*3 

92  CKV!LLBS*L)=RGBM!M2*K*L) 

200  CONTINUE 

201  CONTINUE 

250  CONTINUE 
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C  CHECK  CURRENT  LINE'S  LABELS 

IF  (IOP(3).EQ.l)  WRITE(7*1545)(LABEL(2*K)*K=1*32) 

300  CONTINUE 
900  CONTINUE 

IF  ( IOP (1 )• EQ. 9 )  WRITE (7* 1502 }(CNO(N)*N=i*LLBS) 

IF  (IOP (1  )•  EQ.  9  )  WR ITE( 7*1501 )IMRGE*LLBS 
WRITE! NOFF* 1501) IMRGE*LLBS 
1545  FORMAT!'  LABEL',3212) 

1501  FORMAT(/'  MERGE  ITERATION: '*15* '  END  WITH  LLBS:  '*15/) 
RETURN 
END 
C 

C  - 

C  SORTING  THE  KERNEL  VECTORS 
C 

SUBROUTINE  SORT(NCLRS) 

REAL  TEMP(3) 

COMMON  /BLQCKO/IOPCIO)* NQFl* NOF2* NOF3*NOFC,NOFD* NOFF 
COMMON  /BLOCK2/CKV( 60,3 ),CN 0(60 )*FKV (30,3)* FNO(30) 

DO  30  I=2,NCLRS 
I2=NCLRS+1-I 
DO  20  J=1*I2 

IF  (CKV(J+1*1).G£.CKV(J*1))  GOTO  20 

DO  10  L  =1*3 

TEMP(L)=CKV(J*1*L) 

CKV(J+1*L)=CKV(J*L) 

CKV(J*L)=TEMP(L) 

10  CONTINUE 

20  CONTINUE 

30  CONTINUE 

WRITE! NOFF* 1533 ) 

1533  FGRMAT( /'  SORTING  CANDIDATE  KERNEL  VECTORS'/) 

RETURN 

END 

C 

C  - 

C  TO  FIND  FINAL  KERNEL  VECTORS 

C  LIMIT  TO  10  ITERATIONS 

C 

SUBROUTINE  KERVEC(NC,KK,DIS) 

C  DIST  ARRAY  STORES  THE  TOTAL  DISTANCES  OF  ITERATIONS 

C  C  ARRAY  STORES  NUMBER  OF  PIXELS  FOR  EACH  CLUSTER 

C  D  ARRAY  STORES  TEMPORARY  DISTANCES  TO  CLUSTER  CENTERS 

C  FOR  CURRENT  PIXEL  BEING  PROCESSED 

REAL  DIST (10 )*FIRST(60 ) 

COMMON  /BLOCKO/ IOP(l 0)*NOF1,NOF2*NOF3*NOFC*NOFD*NQFF 
COMMON  /BLOCKl/NOL,NOP*NOL2*  NOP2*NOLH*NOPH 
COMMON  /BLOCK2/CKV (60*3 ),CNO (60), FKV (30*3 ),FNO(30 ) 
COMMON  /BLOCK 3 /IR( 2 56), IG(256 )* IB (2 5 6), RGB (256*3) 
COMMON  /BLOCK5/LAB£L(2* 128), ICOLOR( 256), IBW(256) 

1020  FORMATC  IN  KERVEC*  KK:'*I5) 

C  FINAL  KERNELVECTORS  SAVED  IN  FKV  ARRAY 

C  FNQ  STORE  NO.  OF  PIXELS  IN  EACH  CLUSTER 

C  K-MEANS  ALGORITHM 

C  10-JUM-82  CORRECT  IMPLEMENTATION 
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C  REFERENCE:  TOU  AND  GONZALEZ,**  PATTERN  RECOGNITION 

C  PRINCIPLES",  PP. 94-97 

DO  12  N=1,NC 
FNO(N)-CNOCN) 

CNO(N)=0. 

DO  10  L=l,3 
FKV(N,L  )=CKV(N,L) 

CKV!N,L )-0. 

10  CONTINUE 

12  CONTINUE 

C  NO.  OF  ITERATIONS  LIMIT  TO  10 

DO  500  KK=1, 10 
WRITE! 7,  1020JKK 

WRITE! NOFF,1020 )KK 
DO  15  N=1,NC 
15  FNO!N)=0. 

REWIND  NOFi 
REWIND  NOF2 
REWIND  NOF3 
REWIND  NOFD 

C  FOR  EACH  ITERATION: 

C  REWIND  THE  FEATURES  FILES:  R  ,  G  ,  B  COMPONENTS 

C  REWIND  THE  TEMPORARY  CLASSIFIED  RESULT  FILE,  NOFD 

C  CLASSIFYING  STANDARDS  IN  FRY  ARRAY 

C  AT  THE  SAME  TIME,  COLLECTING  THE  NEW  CENTERS  IN  CKV  ARRAY 

C  I.E.,  UPDATING  THE  KERNEL  YECTORS  BY  CURRENT  CLUSTERING 

C  ICOLOR  ARRAY  STORES  CLASSIFIED  RESULT  OF  CURRENT  LINE 

DO  200  1*1, NOL 
INDX1-I 

READ!N0F1'INDX1)!IR!K),K=1,N0P) 

INDX2=I 

READ!NOF2'INDX2)!IG!K),K=l,NOP) 

INDX3-I 

RSAD!N0F3’1NDX3)!1B!K),K-1,N0P) 

DO  20  J=l,NOP 

RGB! J, 1 )=FLOAT!  IR! J)) 

RGB(J,2)=FLOAT!IG!J)) 

RGB! J,3 )=FLOAT(IB (J )) 

20  CONTINUE 

C  GO  THROUGH  PIXELS  TO  LABEL  THEM  WITH  CLUSTERS 

DO  100  J=1,N0P 
DO  40  N=1,NC 
SUM-0. 

DO  30  L-1,3 

30  SUM=SUM*!RGB!J,L)-FK¥!N,L))**2 

40  FIRST! N )=SGRT!SUM) 

CALL  DISMIN!FIRST,FMIN,NN,NC) 

ICOLOR! J)=NN 
CNG!NN)=CNO!NN)*l. 

C  CURRENT  PIXEL  WAS  FOUND  CLOSER  TO  NN-TH  CLUSTER 

C  THE  FEATURES  SHOULD  BE  INCLUDED  TO  UPDATE  THE  NN-TH 

C  KERNEL  VECTOR 

DO  70  L=l,3 

70  CKV!NN,L)=!CKV!NN,L)*!CNO!NN)-l.)+RGB!J,L))/CNO!NN) 

100  CONTINUE 
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mcxosi 

IF  (10P(2).EQ.l)  MRITE(7,1505)(ICOLORCK),K=1,64) 
HRITECNQFD*INDXD)CICOLORCK),K=l,NOP) 

200  CONTINUE 

C  CURRENT  IN  CKV)  STORE  IN  FKV  TO  BE  USED  TO  CLASSIFY 
C  IN  ROUTINE  DISTAN,  AND  MILL  GIVE  TOTAL  DISTANCE 
DO  260  N*1,NC 
FNO(N)=CNO(N) 

CNO(N)=0. 

DO  250  L=l,3 
FKVCN,L)=CKVCN,L) 

CKVCN,L)=0. 

250  CONTINUE 
260  CONTINUE 

NRITEC NOFF,1503)CCFKVCN,L),L=1,3),N=1,NC) 

1503  FORMAT C /10X, *CH ECK  FKV:  '/ 

130((5X,3E20.8)/)> 

MRITEC  NOFF,1506)CFNOCN),N=1, NC) 

1506  FORMATC*  OF  PIXELS  IN  EACH  CLUSTER:*/ 

1  30C/F12.1 )) 

C  ******** 

CALL  DISTANCNC,DIS) 

C  ******** 

DIST(KK)=DIS 

MRITEC 7, 1504)DISTCKK),KK 
MRITECNOFF,1504)DISTCKK),KK 

1504  FORMATC/*  IN  KERVEC:  DISTCKK)=*,E20. 8, *  KK  =  *,I2/) 
IFCKK.  EQ«  l)GO  TO  500 

RATIO=CDISTCKK)~DISTCKK-l))/DISTCKK-l) 

NRITEC NOFF,1005)RATIQ 
IF  CABS ( RATIO) «LT. 0. 001 )  GOTO  900 
500  CONTINUE 
900  CONTINUE 

1005  FORMATC*  RATIO  IN  KERVEC: *,E20. 8) 

1505  FORMATC*  COLOR: *,6411) 

RETURN 

END 

C  - 

C 

C  USE  CURRENT  KERNEL  VECTORS  TO  CALCULATE  THE  TOTAL 
C  DISTANCE  OF  THE  IMAGE 

C  NC:  NUMBER  OF  CLUSTERS 

C  DISTOT:  C  RESULT  )  TOTAL  DISTANCE 
C 

SUBROUTINE  DISTANCNC,DISTOT) 

COMMON  /BLOCKO/ IOP CIO), NOF1 , NOF2, N0F3, NOFC, NOFD, NOFF 
COMMON  /BLOCK 1/NOL, NOP, NOL2, NOP2,NOLH, NOPH 
COMMON  /BLOCK2/CKVC60,3),CNOC60),FKVC30,3),FNOC30) 
COMMON  /BLOCK3/IRC256), IGC256),IBC256), RGB C 256, 3) 
COMMON  /BL0CK5/LABELC2,128),IC0L0RC256),IBNC256) 
DISTQf =0. 

REHIND  NOF1 
REMIND  NOF2 
REMIND  NOF3 
REMIND  NOFD 
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DO  200  I=l,NOL 
INDXD~I 

READ!NOFD'INDXD)!ICOLOR!K),K=1,NQP) 

INDXl'I 

READINOFl'INDXl )!IR!K),K=l,NOP) 

INDX2-I 

READ!NQF2'INDX2 )!IG!K),K>1, NOP) 

INDX3=I 

READ!N0F3'INDX3)!IB!K),K=1,N0P) 

C  STORE  FEATURE  VECTOR  IN  WORKING  VARIABLE  X 

DO  10  K=l,NOP 
RGBCK/l  )=PLOAT(IR(K)) 

RGB(K/2 )=FLOAT(IG(K )) 

RGB!K,3  )=FLOAT(IB(K)) 

10  CONTINUE 

DO  100  J=l,NOP 
NCLSR= I COLOR ( J ) 

SUH=0. 

DO  30  L-l/3 

SUM=SUM+(RGB! J,L)-FKV! NCLSR/L ))**2 
30  CONTINUE 

D I STOTsDISTOT-frSQRT ( SUM ) 

100  CONTINUE 

200  CONTINUE 

NR ITE( NOFF,1501 JDISTOT 

1501  FORMAT!/*  IN  D1STAN:  DISTOT  *  *,E20.8/> 

RETURN 
END 

C 
C 

C  USE  THE  FINAL  KERNEL  VECTORS  TO  CLASSIFY  THE  IMAGE 

C  INTO  CLUSTERS  (  SUBREGIONS  ) 

C 

SUBROUTINE  CLASS(NC) 

REAL  FIRST(60) 

COMMON  /BLOCKO/ IOP! 10),NOF1,NOF2,NOF3,NOFC,NOFD,NOFF 
COMMON  /BLOCK 1 / NOL, NOP , NOL2 / NOP2 , NOL H, NOPH 
COMMON  /BLOCK2/CKV!60,3),CNO!60),FKV!30,3),FNO!30) 
COMMON  /BLOCK3/IR!256),IG!256),IB!256),RGB!256,3) 
COMMON  /BL0CK5/LABEL!2,128),IC0LQR!256),1BW!256) 

WRITE! NOFF, 1505) 

1505  FORMAT! /*  IN  CLASS:  V) 

WRITE! 7, 1501 )!!FKV!N,L),L=l,3),Nsl,NC) 

WRITE! NOFF,1501)!!FKV!N,L),L=1,3),N=1,NC) 

1501  FORMAT! 10 X/'  FINAL  KERNEL  VECTERS  ! */!30! 5X,3E20.8/ ))) 
1509  FORMAT! IX, 6411) 

C  USE  FINAL  KERNAL  VECTORS  TO  CLASSIFY  THE  PICTURE 
C  NOFC:  INTEGER  OUTPUT)  NOFD!  REAL  OUTPUT 
REWIND  NOF1 
REWIND  NOF2 
REWIND  NOF3 
REWIND  NOFC 
REWIND  NOFD 

C  RANGE  OF  FINAL  RESULT  0..180 
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FACT=1 80. /FLOAT (NC) 

DO  200  Isl,NOL 
INDX1-I 
INDX2-I 
INDX3-I 

READCNOFl'INDXl )(IR(K)#K=l#NOP) 
READ(NOF2'INDX2)(IG(K)#K=l#NOP) 

READ(NOF3*INDX3 )(IB(K)#K=l#NOP) 

DO  10  J=l#NOP 
RGB(J,l)=FLOAT(IR(J)) 

RGB( J/2 )=FLOAT( IG( J )) 

RGB( J#  3 )=FLOAT( IB ( J ) ) 

10  CONTINUE 

DO  100  J-l,NOP 
DO  40  N=1#NC 
SUM=0. 

DO  30  L=L,3 

SUM=SUM+( RGB( J#L)-FKV(N#L ))**2 
30  CONTINUE 

FIRST(N)~SGRT(SUM) 

40  CONTINUE 

C  DISTANCES  TO  FINAL  KERNEL  VECTOS  OF  NC  CLUSTERS 

C  FROM  CURRENT  PIXEL  ARE  STORED  IN  DIS  ARRAY# 

C  CALLING  SUBROUTINE  DISM1N  TO  FIND  TO  HHICH  CLUSTER 

C  THE  CURRENT  PIXEL  IS  CLOSER  (  MINIMUM  DISTANCE  ). 

C  THE  RESULT  IS  KMIN-TH  CLUSTER 

CALL  DISMIN(FIRST#SMIN#NMIN#NC) 

C  BLACK  AND  WHITE  DISPLAY  PURPOSE:  NEEDS  TO  MULTIPLY  A  FACTOR 
C  TO  BE  IN  REASONABLE  CRAY  LEVEL  RANGE 

IBW(J) =INT(FLOAT(NMIN)*FACT) 

C  COLOR  DISPLAY  PURPOSE:  AN  INTEGER  IN  RANGE  1  TO  7 
IC0L0R(J)-NM1N 
100  CONTINUE 

INDXD-I 

NR  ITE(  NOf  D'INDXDH  IBW(K)#K=1#  NOP) 

INDXC-I 

HR ITE( NOFC'INDXC)(ICOLOR(K)#K=l#NOP) 

IF  (I.LE.64)  WRITE(7#1509)(ICOLOR(K)#K=1#64) 

IF  (I.LE.64)  WRITE(N0FF#15Q9) (ICQL0R(K)#K=1#64) 

200  CONTINUE 

RETURN 
END 
C 

C 

SUBROUTINE  DISMIN(DARRY#DATMIN#NOMIN#NCLSTR) 

C  PASS  DARRY  ARRAY  MITH  NCLSTR  ELEMENTS  MEANINGFUL 
C  SEARCH  FOR  THE  MINIMUM  ELEMENT#  NQMIN-TH  ELEMENT# 

C  HITH  VALUE  DATMIN;  PASS  BACK  D ATM IN  AND  NOMIN  BACK 
C  CALLING  ROUTINE 

REAL  DARRYC60) 

DATMIN sDARRY(l) 

NOMIN-1 

C  ASSUME  THE  FIRST  ELEMENT  IS  THE  MINIMUM 

C  THEN  GO  THROUGH  THE  REST  OF  THE  ARRAY  TO  FIND  ANY  SMALLER 


■ 
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IF  (NCLSTft.EQ.l)  GOTO  900 
DO  100  I=2,NCLSTR 
IF  (DARRy(I).GE.DATMIN)  GOTO  100 
DATMINsDARRYII) 

NOMIN=I 

100  CONTINUE 

900  CONTINUE 
RETURN 
END 


